library(plm)
library(multiwayvcov)
library(lmtest)
library(stargazer)
library(ggplot2)
library(clusterSEs)
library(dplyr)

# setwd("~/Dropbox/Replication_JPT")
data <- read.csv("jpt_south_data.csv", stringsAsFactors = FALSE, header = T)

'%!in%' = Negate('%in%')

# ---------------------
#       Figure 1
# ---------------------

# Figure 1 (a) Taxes per White Capita
fig1a_data <- data %>% 
  filter(!is.na(state_taxes_PWC_5y)) %>%
  group_by(malapportion, year) %>% 
  summarise(mean_group = mean(state_taxes_PWC_5y, na.rm = TRUE))
fig1a <- ggplot(fig1a_data, aes(x = year, y = mean_group)) +
  geom_line(aes(x = year, y = mean_group, color = as.factor(malapportion))) +
  scale_color_discrete(name  ="", breaks=c("0", "1"), labels=c("Non-malapportioned States", "Malapportioned States")) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="bottom") +
  xlim(c(1844, 1860)) + labs(x = "Year", y = "") 
ggsave("fig1a.png", width = 5, height = 5, dpi = 300)

# Figure 1 (b) Taxes/Income
fig1b_data <- data %>% 
  filter(!is.na(st_taxes_per_income_5y)) %>%
  group_by(malapportion, year) %>% 
  summarise(mean_group = mean(st_taxes_per_income_5y))
fig1b <- ggplot(fig1b_data, aes(x = year, y = mean_group)) +
  geom_line(aes(x = year, y = mean_group, color = as.factor(malapportion))) +
  scale_color_discrete(name  ="", breaks=c("0", "1"), labels=c("Non-malapportioned States", "Malapportioned States")) +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="bottom") +
  xlim(c(1844, 1860)) + labs(x = "Year", y = "")
ggsave("fig1b.png", width = 5, height = 5, dpi = 300)

# Figure 1 (c) Ad Valorem Land Tax
fig1c_data <- data %>% 
  filter(!is.na(tax_rate_land_5y)) %>%
  group_by(malapportion, year) %>% 
  summarise(mean_group = mean(tax_rate_land_5y, na.rm = TRUE))
fig1c <- ggplot(fig1c_data, aes(x = year, y = mean_group)) +
  geom_line(aes(x = year, y = mean_group*100, color = as.factor(malapportion))) +
  scale_color_discrete(name  ="", breaks=c("0", "1"), labels=c("Non-malapportioned States", "Malapportioned States")) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="bottom")  +
  xlim(c(1844, 1860)) + ylim(c(12, 29)) + labs(x = "Year", y = "")
ggsave("fig1c.png", width = 5, height = 5, dpi = 300)

# Figure 1 (d) Tax Rate on Slaves
fig1d_data <- data %>% 
  filter(!is.na(slaves_tax_rate_5y)) %>%
  group_by(malapportion, year) %>% 
  summarise(mean_group = mean(slaves_tax_rate_5y, na.rm = TRUE))
fig1d <- ggplot(fig1d_data, aes(x = year, y = mean_group)) +
  geom_line(aes(x = year, y = mean_group*100, color = as.factor(malapportion))) +
  scale_color_discrete(name  ="", breaks=c("0", "1"), labels=c("Non-malapportioned States", "Malapportioned States")) +
  theme_bw() +
  scale_x_continuous(breaks=seq(1845,1860,3)) + ylim(c(9, 15)) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="bottom")  +
  labs(x = "Year", y = "")
ggsave("fig1d.png", width = 5, height = 5, dpi = 300)

# Figure 1 (e) Poll Tax Rates
fig1e_data <- data %>% 
  filter(!is.na(poll_tax_per_wh_male_5y)) %>%
  group_by(malapportion, year) %>% 
  summarise(mean_group = mean(poll_tax_per_wh_male_5y, na.rm = TRUE))
fig1e <- ggplot(fig1e_data, aes(x = year, y = mean_group)) +
  geom_line(aes(x = year, y = mean_group, color = as.factor(malapportion))) +
  scale_color_discrete(name  ="", breaks=c("0", "1"), labels=c("Non-malapportioned States", "Malapportioned States")) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position="bottom")  +
  xlim(c(1844, 1860)) + ylim(c(0, .6)) + labs(x = "Year", y = "")
ggsave("fig1e.png", width = 5, height = 5, dpi = 300)


# -----------------------------------------------------------------------------
#  Table 1: Two-way Fixed Effects: Commodity Price Index and Taxation Outcomes
# -----------------------------------------------------------------------------

pdata <- pdata.frame(data[which(data$year >= 1840),], index = c("statefips", "year"))
base.form <- formula(. ~ log(index.rel_5y)*malapportion + log(tot_pop+1) + log(urban_pop+1) + log(income_linear))

tw.inc <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.inc <- coeftest(tw.inc, vcovHC(tw.inc, type = 'HC0', cluster = 'group'))[,2]
coef.tw.inc <- coeftest(tw.inc, vcovHC(tw.inc, type = 'HC0', cluster = 'group'))

tw.pwc <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.pwc <- coeftest(tw.pwc, vcovHC(tw.pwc, type = 'HC0', cluster = 'group'))[,2]
coef.tw.pwc <- coeftest(tw.pwc, vcovHC(tw.pwc, type = 'HC0', cluster = 'group'))

tw.slv <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.slv <- coeftest(tw.slv, vcovHC(tw.slv, type = 'HC0', cluster = 'group'))[,2]
coef.tw.slv <- coeftest(tw.slv, vcovHC(tw.slv, type = 'HC0', cluster = 'group'))

tw.land <- plm(update(base.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.land <- coeftest(tw.land, vcovHC(tw.land, type = 'HC0', cluster = 'group'))[,2]
coef.tw.land <- coeftest(tw.land, vcovHC(tw.land, type = 'HC0', cluster = 'group'))

tw.poll <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.poll <- coeftest(tw.poll, vcovHC(tw.poll, type = 'HC0', cluster = 'group'))[,2]
coef.tw.poll <- coeftest(tw.poll, vcovHC(tw.poll, type = 'HC0', cluster = 'group'))

stargazer(tw.inc, tw.pwc, tw.slv, tw.land, tw.poll,
          se = list(se.tw.inc, se.tw.pwc, se.tw.slv, se.tw.land, se.tw.poll),
          keep = "index.rel_5y",
          title = "Two-way Fixed Effects with Commodity Price Index",
          covariate.labels = c("Commodity Index (log)", "Commodity Index (log) $times$ Malapportionment"),
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"), 
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"), 
                           c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Number of States","14", "14", "14", "14", "14")) )

set.seed(65398)

# Wild-cluster bootstrap p-values
tw.inc.wild <- cluster.wild.plm(mod=tw.inc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.inc.wild)[5]
tw.pwc.wild <- cluster.wild.plm(mod=tw.pwc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.pwc.wild)[5]
tw.slv.wild <- cluster.wild.plm(mod=tw.slv, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.slv.wild)[5]
tw.land.wild <- cluster.wild.plm(mod=tw.land, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.land.wild)[5]
tw.poll.wild <- cluster.wild.plm(mod=tw.poll, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.poll.wild)[5]

# ------------------------------------------------
#  Figures A1 and A2: see "jpt_jop_maps_2022.R"
# ------------------------------------------------

# ---------------------
#       Figure A3
# ---------------------

##  Figure A3 Balance Test for Pre-Treatment Covariates, 1840
df <- data %>%
  filter(year == 1840) %>%
  select(stateabbrev, malapportion, sugarsuit, tobaccosuit, cottonsuit, tot_pop, pwc.public.officials1850, 
         urban.pc, slave.pc, income_linear, rivers_area, cotton.prod1840, qcd)

t.1 <- t.test(slave.pc*100 ~ malapportion, data = df)$estimate[1]
t.2 <- t.test(urban.pc*100 ~ malapportion, data = df)$estimate[1] 
t.3 <- t.test(income_linear/1000000 ~ malapportion, data = df)$estimate[1] 
t.4 <- t.test(rivers_area*100 ~ malapportion, data = df)$estimate[1]
t.5 <- t.test(cottonsuit ~ malapportion, data = df)$estimate[1]
t.6 <- t.test(pwc.public.officials1850*100 ~ malapportion, data = df)$estimate[1]
t.7 <- t.test(sugarsuit ~ malapportion, data = df)$estimate[1]
t.8 <- t.test(tobaccosuit ~ malapportion, data = df)$estimate[1]
t.9 <- t.test(cotton.prod1840/1000000 ~ malapportion, data = df)$estimate[1]
t.10 <- t.test(qcd ~ malapportion, data = df)$estimate[1]

g.1 <- t.test(slave.pc*100 ~ malapportion, data = df)$estimate[2]
g.2 <- t.test(urban.pc*100 ~ malapportion, data = df)$estimate[2] 
g.3 <- t.test(income_linear/1000000 ~ malapportion, data = df)$estimate[2]
g.4 <- t.test(rivers_area*100 ~ malapportion, data = df)$estimate[2]
g.5 <- t.test(cottonsuit ~ malapportion, data = df)$estimate[2]
g.6 <- t.test(pwc.public.officials1850*100 ~ malapportion, data = df)$estimate[2]
g.7 <- t.test(sugarsuit ~ malapportion, data = df)$estimate[2]
g.8 <- t.test(tobaccosuit ~ malapportion, data = df)$estimate[2]
g.9 <- t.test(cotton.prod1840/1000000  ~ malapportion, data = df)$estimate[2]
g.10 <- t.test(qcd ~ malapportion, data = df)$estimate[2]

p.1 <- t.test(slave.pc*100 ~ malapportion, data = df)$p.value
p.2 <- t.test(urban.pc*100 ~ malapportion, data = df)$p.value
p.3 <- t.test(income_linear/1000000 ~ malapportion, data = df)$p.value
p.4 <- t.test(rivers_area*100 ~ malapportion, data = df)$p.value
p.5 <- t.test(cottonsuit ~ malapportion, data = df)$p.value
p.6 <- t.test(pwc.public.officials1850*100 ~ malapportion, data = df)$p.value
p.7 <- t.test(sugarsuit ~ malapportion, data = df)$p.value
p.8 <- t.test(tobaccosuit ~ malapportion, data = df)$p.value
p.9 <- t.test(cotton.prod1840/1000000 ~ malapportion, data = df)$p.value
p.10 <- t.test(qcd ~ malapportion, data = df)$p.value

mean.g1 <- as.data.frame(rbind(t.1, t.2, t.3, t.4, t.5, t.6, t.7, t.8, t.9, t.10)) 
mean.g2 <- as.data.frame(rbind(g.1, g.2, g.3, g.4, g.5, g.6, g.7, g.8, g.9, g.10)) 
p.values <- as.data.frame(rbind(p.1, p.2, p.3, p.4, p.5, p.6, p.7, p.8, p.9, p.10)) 
diffs <- cbind(mean.g1, mean.g2, p.values)
names(diffs)[3] <- "p.value"
names(diffs)[1] <- "mean_0"
names(diffs)[2] <- "mean_1"
diffs$model <- c("Enslaved Population (%)", "Urban Population (%)", "State Income ($ millions)", 
                 "Density of Navigable Rivers", "Cotton Suitability Index", 
                 "Public Officials 1850 (% white pop)", "Sugar Suitability Index", "Tobacco Suitability Index",
                 "Cotton Production, 1840 (million pounds)", "Dispersion in Cotton Suitability")

diffs$mean_0 <- format(round(diffs$mean_0, 2), nsmall = 2)
diffs$mean_1 <- format(round(diffs$mean_1, 2), nsmall = 2)
diffs <- diffs[order(diffs$p.value), ]  
diffs$model <- factor(diffs$model, levels = diffs$model)

# Balance Plot Code by Rocio Titiunik 
# https://titiunik.mycpanel.princeton.edu/software/balance-plots/graph.pval.public.R
plot.pval <- function(results, title=NULL, legend,legendx=0.15,legendy=2.2, textsize=0.9, parcex=0.8, at1=-0.35, at2=-0.15, at3=-0.9,xlim1=-0.85) {
  # set values of different parameters
  xlim = c(xlim1,1); pchset = c(21,24,22,23); pchcolset = c("blue","yellow","red","darkgreen")
  # set margins and letter size
  par(cex=parcex, mai = c(0.5, 0.35, 1.1, 0.35))
  # set number of rows 
  ny = nrow(results)
  # create the empty figure
  if(!is.null(title))  plot(x=NULL,axes=F, xlim=xlim, ylim=c(1,ny),xlab="",ylab="", main=title)
  if(is.null(title))   plot(x=NULL,axes=F, xlim=xlim, ylim=c(1,ny),xlab="",ylab="")
  # add the 0, 0.05 and 0.1 vertical lines
  abline(v=c(0,0.05,0.1),lty=c(1,4,4), lwd=c(1,2,2))
  axis(side=1,at=c(0,0.05,0.1,1),tick=TRUE, las=2, cex.axis=0.7)
  # add labels on top of the three areas of the graph
  axis(side=3,at=at1,labels="Mean\nMS",tick=FALSE, padj=0.5,cex.axis=textsize)
  axis(side=3,at=at2,labels="Mean\nNMS",tick=FALSE, padj=0.5,cex.axis=textsize)
  axis(side=3,at=0.5,labels="P-values",tick=FALSE, padj=0.5,cex.axis=textsize)
  # Fill the figure with the information which is inside the 'results' matrix
  # First, add the p-values as points
  for(i in 4:ncol(results)) points(results[,i],ny:1, pch = pchset[i-4+1], col = pchcolset[i-4+1], bg = pchcolset[i-4+1])
  # Second, add each variable name and the means for treated and control
  for(i in 1:ny) {
    text(at3,ny-i+1,results[i,1],adj = 0,cex=textsize) # variable name
    text(at1,ny-i+1,results[i,2], cex=textsize)        # treatment mean
    text(at2,ny-i+1,results[i,3], cex=textsize)        # control mean
  }
  # Add dotted horizontal lines every two variables to make it prettier
  for(i in seq(2,by=2,length.out=floor((ny-1)/2))) abline(h = i+0.5, lty = 3)
  # Add legend
  if(legend) legend(x=legendx, y=legendy, c("Diff-means", "Diff-medians", "D-stat", "Rank-sum"), pch=pchset, pt.bg = pchcolset, cex=0.8)
}

results <- diffs[, c("model", "mean_1", "mean_0", "p.value")]
# png(file = "balance.png", width=10, height=5, units="in", res = 100)
par(mfrow=c(1,1))
plot.pval(results=results, title="", legend=FALSE,legendx=0.15,legendy=2.2, textsize=0.9, parcex=0.8, at1=-0.35, at2=-0.15, at3=-0.9,xlim1=-0.85)
# dev.off()

# ---------------------
#       Figure A4
# ---------------------

# Figure A4: Marginal Effects of Commodity Price Index (log)
dat <- data[,c("st_taxes_per_income_5y", "state_taxes_PWC_5y", "tax_rate_land_5y", "log.slaves_tax_rate_5y", "poll_tax_per_wh_male_5y")]

# Column 1: TWFE w/o controls
for (i in 1:5){
  outcome = colnames(dat)[i]
  mod1 <- plm(get(outcome) ~ log.index.rel_5y:malapportion + log.index.rel_5y, data = data, index=c("stateabbrev", "year"), model="within", effect = "twoways")
  beta.hat <- coef(mod1) 
  cov <- vcovHC(mod1, type = 'HC0', cluster = 'group')
  z0 <- c(0, 1)
  dy.dx <- beta.hat["log.index.rel_5y"] + beta.hat["log.index.rel_5y:malapportion"]*z0
  se.dy.dx <- sqrt(cov["log.index.rel_5y", "log.index.rel_5y"] + z0^2*cov["log.index.rel_5y:malapportion", "log.index.rel_5y:malapportion"] + 2*z0*cov["log.index.rel_5y", "log.index.rel_5y:malapportion"])
  upr <- dy.dx + 1.96*se.dy.dx
  lwr <- dy.dx - 1.96*se.dy.dx
  temp_plot = ggplot(data=NULL,aes(x=z0, y=dy.dx)) + 
    labs(x="Malapportionment", y="Marginal Effect",
         title="") +
    geom_pointrange(aes(ymin=lwr, ymax=upr)) +
    geom_hline(yintercept=0, linetype = "dashed") +
    theme_classic() + scale_x_discrete(limits=c(0, 1)) +
    theme(axis.title.x = element_text(size=20),
          axis.title.y = element_text(size=20) )
  ggsave(temp_plot, file=paste0("fe_mod1_", i,".png"), bg = "white")
}
# Figures: fe_mod1_1 = "Taxes/Income"; fe_mod1_2 = "Taxes PWC"; fe_mod1_3 = "Tax Rate on Land"; fe_mod1_4 = "Tax Rate on Slaves (log)"; fe_mod1_5 = "Poll Tax Rate"

# Column 2: TWFE w/ controls
for (i in 1:5){
  outcome = colnames(dat)[i]
  mod2 <- plm(get(outcome) ~ log.index.rel_5y:malapportion + log.index.rel_5y + log(tot_pop+1) + log(urban_pop+1) + log(income_linear+1), data = data, index=c("stateabbrev", "year"), model="within", effect = "twoways")
  beta.hat <- coef(mod2) 
  cov <- vcovHC(mod2, type = 'HC0', cluster = 'group')
  z0 <- c(0, 1)
  dy.dx <- beta.hat["log.index.rel_5y"] + beta.hat["log.index.rel_5y:malapportion"]*z0
  se.dy.dx <- sqrt(cov["log.index.rel_5y", "log.index.rel_5y"] 
                   + z0^2*cov["log.index.rel_5y:malapportion", "log.index.rel_5y:malapportion"] 
                   + 2*z0*cov["log.index.rel_5y", "log.index.rel_5y:malapportion"])
  upr <- dy.dx + 1.96*se.dy.dx 
  lwr <- dy.dx - 1.96*se.dy.dx
  temp_plot = ggplot(data=NULL,aes(x=z0, y=dy.dx)) + 
    labs(x="Malapportionment", y="Marginal Effect",
         title="") +
    geom_pointrange(aes(ymin=lwr, ymax=upr)) +
    geom_hline(yintercept=0, linetype = "dashed") +
    theme_classic() + scale_x_discrete(limits=c(0, 1)) +
    theme(axis.title.x = element_text(size=20),
          axis.title.y = element_text(size=20) )
  ggsave(temp_plot, file=paste0("fe_mod2_", i,".png"), bg = "white")
}
# Figures: fe_mod2_1 = "Taxes/Income"; fe_mod2_2 = "Taxes PWC"; fe_mod2_3 = "Tax Rate on Land"; fe_mod2_4 = "Tax Rate on Slaves (log)"; fe_mod2_5 = "Poll Tax Rate"


# Column 3: First-Difference w/o controls
for (i in 1:5){
  outcome = colnames(dat)[i]
  fd1 <- plm(get(outcome) ~ log.index.rel_5y:malapportion + log.index.rel_5y, data = data, index=c("stateabbrev", "year"), model="fd")
  beta.hat <- coef(fd1) 
  cov <- vcov(fd1)
  z0 <- c(0, 1)
  dy.dx <- beta.hat["log.index.rel_5y"] + beta.hat["log.index.rel_5y:malapportion"]*z0
  se.dy.dx <- sqrt(cov["log.index.rel_5y", "log.index.rel_5y"] + z0^2*cov["log.index.rel_5y:malapportion", "log.index.rel_5y:malapportion"] + 2*z0*cov["log.index.rel_5y", "log.index.rel_5y:malapportion"])
  upr <- dy.dx + 1.96*se.dy.dx
  lwr <- dy.dx - 1.96*se.dy.dx
  temp_plot = ggplot(data=NULL,aes(x=z0, y=dy.dx)) + 
    labs(x="Malapportionment", y="Marginal Effect",
         title="") +
    geom_pointrange(aes(ymin=lwr, ymax=upr)) +
    geom_hline(yintercept=0, linetype = "dashed") +
    theme_classic() + scale_x_discrete(limits=c(0, 1)) +
    theme(axis.title.x = element_text(size=20),
          axis.title.y = element_text(size=20) )
  ggsave(temp_plot, file=paste0("fd_mod1_", i,".png"), bg = "white")
}
# Figures: fd_mod1_1 = "Taxes/Income"; fd_mod1_2 = "Taxes PWC"; fd_mod1_3 = "Tax Rate on Land"; fd_mod1_4 = "Tax Rate on Slaves (log)"; fd_mod1_5 = "Poll Tax Rate"


# Column 4: First-Difference w/ controls
for (i in 1:5){
  outcome = colnames(dat)[i]
  fd2 <- plm(get(outcome) ~ log.index.rel_5y:malapportion + log.index.rel_5y + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = data, index=c("stateabbrev", "year"), model="fd")
  beta.hat <- coef(fd2) 
  cov <- vcov(fd2)
  z0 <- c(0, 1)
  dy.dx <- beta.hat["log.index.rel_5y"] + beta.hat["log.index.rel_5y:malapportion"]*z0
  se.dy.dx <- sqrt(cov["log.index.rel_5y", "log.index.rel_5y"] + z0^2*cov["log.index.rel_5y:malapportion", "log.index.rel_5y:malapportion"] + 2*z0*cov["log.index.rel_5y", "log.index.rel_5y:malapportion"])
  upr <- dy.dx + 1.96*se.dy.dx
  lwr <- dy.dx - 1.96*se.dy.dx
  temp_plot = ggplot(data=NULL,aes(x=z0, y=dy.dx)) + 
    labs(x="Malapportionment", y="Marginal Effect",
         title="") +
    geom_pointrange(aes(ymin=lwr, ymax=upr)) +
    geom_hline(yintercept=0, linetype = "dashed") +
    theme_classic() + scale_x_discrete(limits=c(0, 1)) +
    theme(axis.title.x = element_text(size=20),
          axis.title.y = element_text(size=20) )
  ggsave(temp_plot, file=paste0("fd_mod2_", i,".png"), bg = "white")
}
# Figures: fd_mod2_1 = "Taxes/Income"; fd_mod2_2 = "Taxes PWC"; fd_mod2_3 = "Tax Rate on Land"; fd_mod2_4 = "Tax Rate on Slaves (log)"; fd_mod2_5 = "Poll Tax Rate"


# -----------------------------------------------------
#     Table A2: Cotton Prices and Taxation Outcomes
# -----------------------------------------------------
cotton.form <- formula(. ~ log(cotton.newOrleans_5y)*malapportion + log(tot_pop+1) + log(urban_pop+1) + log(income_linear))

tw.inc.cot <- plm(update(cotton.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.inc.cot <- coeftest(tw.inc.cot, vcovHC(tw.inc.cot, type = 'HC0', cluster = 'group'))[,2]

tw.pwc.cot <- plm(update(cotton.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.pwc.cot <- coeftest(tw.pwc.cot, vcovHC(tw.pwc.cot, type = 'HC0', cluster = 'group'))[,2]

tw.slv.cot <- plm(update(cotton.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.slv.cot <- coeftest(tw.slv.cot, vcovHC(tw.slv.cot, type = 'HC0', cluster = 'group'))[,2]

tw.land.cot <- plm(update(cotton.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.land.cot <- coeftest(tw.land.cot, vcovHC(tw.land.cot, type = 'HC0', cluster = 'group'))[,2]

tw.poll.cot <- plm(update(cotton.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.tw.poll.cot <- coeftest(tw.poll.cot, vcovHC(tw.poll.cot, type = 'HC0', cluster = 'group'))[,2]

stargazer(tw.inc.cot, tw.pwc.cot, tw.slv.cot, tw.land.cot, tw.poll.cot,
          se = list(se.tw.inc.cot, se.tw.pwc.cot, se.tw.slv.cot, se.tw.land.cot, se.tw.poll.cot), 
          keep = "cotton.newOrleans_5y",
          covariate.labels = c("Cotton Prices $times$ Malapportionment"),
          title = "Cotton Prices and Taxation Outcomes",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"), 
                           c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Number of States","14", "14", "14", "14", "14")) )

set.seed(63598)
# Wild-cluster bootstrap p-values
tw.inc.wild.cot <- cluster.wild.plm(mod=tw.inc.cot, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.inc.wild.cot)[4]
tw.pwc.wild.cot <- cluster.wild.plm(mod=tw.pwc.cot, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.pwc.wild.cot)[4]
tw.slv.wild.cot <- cluster.wild.plm(mod=tw.slv.cot, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.slv.wild.cot)[4]
tw.land.wild.cot <- cluster.wild.plm(mod=tw.land.cot, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.land.wild.cot)[4]
tw.poll.wild.cot <- cluster.wild.plm(mod=tw.poll.cot, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(tw.poll.wild.cot)[4]


# --------------------------------------------------------------------------------------
#    Table A3: Time-Invariant Controls: Commodity Price Index and Taxation Outcomes
# --------------------------------------------------------------------------------------
inv.form <- formula(. ~ log(index.rel_5y)*malapportion + log(tot_pop+1) + log(urban_pop+1) + log(income_linear) + rivers_length*as.factor(year) + coastal*as.factor(year) + area*as.factor(year) + cottonsuit*as.factor(year))

tw.inc.i <- plm(update(inv.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.inc.i <- coeftest(tw.inc.i, vcovHC(tw.inc.i, type = 'HC0', cluster = 'group'))[,2]

tw.pwc.i <- plm(update(inv.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.pwc.i <- coeftest(tw.pwc.i, vcovHC(tw.pwc.i, type = 'HC0', cluster = 'group'))[,2]

tw.slv.i <- plm(update(inv.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.slv.i <- coeftest(tw.slv.i, vcovHC(tw.slv.i, type = 'HC0', cluster = 'group'))[,2]

tw.land.i <- plm(update(inv.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.land.i <- coeftest(tw.land.i, vcovHC(tw.land.i, type = 'HC0', cluster = 'group'))[,2]

tw.poll.i <- plm(update(inv.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.poll.i <- coeftest(tw.poll.i, vcovHC(tw.poll.i, type = 'HC0', cluster = 'group'))[,2]

stargazer(tw.inc.i, tw.pwc.i, tw.slv.i, tw.land.i, tw.poll.i,
          se = list(se.inc.i, se.pwc.i, se.slv.i, se.land.i, se.poll.i), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index (log) $times$ Malapportionment"),
          title = "Time-Invariant Controls: Commodity Price Index and Taxation Outcomes",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"), 
                           c("State FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Time-invariant controls $times$ Year FE", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Number of States","14", "14", "14", "14", "14")) )

# Wild-cluster bootstrap p-values
p.inc.wild <- cluster.wild.plm(mod=tw.inc.i, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
p.pwc.wild <- cluster.wild.plm(mod=tw.pwc.i, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
p.slv.wild <- cluster.wild.plm(mod=tw.slv.i, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
p.land.wild <- cluster.wild.plm(mod=tw.land.i, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
p.poll.wild <- cluster.wild.plm(mod=tw.poll.i, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]

# ---------------------------------------------------------------
#     Figure A5: Parallel Trends in State Taxes (% Income)
# ---------------------------------------------------------------
library(npregfast)

# Taxes/Income in the NMS
fit0 <- frfast(st_taxes_per_income ~ year, data = data[which(data$malapportion == 0),], 
               nboot = 1000, na.action = "na.omit", 
               model = "np", smooth = "kernel", kernel = "epanech",)
dataplot0 <- as.data.frame(fit0$pl)[1]
names(dataplot0)[1] <- "pl"
pu <- as.data.frame(fit0$pu)[1]
names(pu) <- "pu"
p <-  as.data.frame(fit0$p)[1]
names(p) <- "p"
x <- as.data.frame(fit0$x)
names(x) <- "x"
dataplot0 <- cbind(dataplot0, pu, p, x)
dataplot0$malapportion <- 0

# Taxes/Income in the MS
fit1 <- frfast(st_taxes_per_income ~ year, data = data[which(data$malapportion == 1),], 
               nboot = 1000, na.action = "na.omit", smooth = "kernel", kernel = "epanech",)
dataplot1 <- as.data.frame(fit1$pl)[1]
names(dataplot1)[1] <- "pl"
pu <- as.data.frame(fit1$pu)[1]
names(pu) <- "pu"
p <-  as.data.frame(fit1$p)[1]
names(p) <- "p"
x <- as.data.frame(fit1$x)
names(x) <- "x"
dataplot1 <- cbind(dataplot1, pu, p, x)
dataplot1$malapportion <- 1

dataplot <- rbind(dataplot0, dataplot1)

# Fit lines w/ confidence intervals
figA5 <- ggplot(dataplot, aes(x = x, y = p)) +
  geom_line(aes(y = p, color = as.factor(malapportion))) +
  geom_line(aes(y = pl, color = as.factor(malapportion)), linetype = 2) +
  geom_line(aes(y = pu, color = as.factor(malapportion)), linetype = 2) +
  scale_color_discrete(name  ="States", breaks=c("0", "1"), labels=c("Non-malapportioned", "Malapportioned")) +
  theme_classic() + theme(legend.position="bottom") +
  labs(x = "Year", y = "State Taxes (% Income)")
figA5
#ggsave("figA5.png", width = 5, height = 5, dpi = 300)

# ---------------------------------------------------------------
#     Figure A6: Parallel Trends in Taxes Per White Capita
# ---------------------------------------------------------------

# Taxes PWC in the NMS
fit0 <- frfast(state_taxes_PWC ~ year, data = data[which(data$malapportion == 0),], 
               nboot = 1000, na.action = "na.omit", 
               model = "np", smooth = "kernel", kernel = "epanech",)
dataplot0 <- as.data.frame(fit0$pl)[1]
names(dataplot0)[1] <- "pl"
pu <- as.data.frame(fit0$pu)[1]
names(pu) <- "pu"
p <-  as.data.frame(fit0$p)[1]
names(p) <- "p"
x <- as.data.frame(fit0$x)
names(x) <- "x"
dataplot0 <- cbind(dataplot0, pu, p, x)
dataplot0$malapportion <- 0

# Taxes PWC in the MS
fit1 <- frfast(state_taxes_PWC ~ year, data = data[which(data$malapportion == 1),], 
               nboot = 1000, na.action = "na.omit", smooth = "kernel", kernel = "epanech",)
dataplot1 <- as.data.frame(fit1$pl)[1]
names(dataplot1)[1] <- "pl"
pu <- as.data.frame(fit1$pu)[1]
names(pu) <- "pu"
p <-  as.data.frame(fit1$p)[1]
names(p) <- "p"
x <- as.data.frame(fit1$x)
names(x) <- "x"
dataplot1 <- cbind(dataplot1, pu, p, x)
dataplot1$malapportion <- 1

dataplot <- rbind(dataplot0, dataplot1)

# Fit lines w/ confidence intervals
figA6 <- ggplot(dataplot, aes(x = x, y = p)) +
  geom_line(aes(y = p, color = as.factor(malapportion))) +
  geom_line(aes(y = pl, color = as.factor(malapportion)), linetype = 2) +
  geom_line(aes(y = pu, color = as.factor(malapportion)), linetype = 2) +
  scale_color_discrete(name  ="States", breaks=c("0", "1"), labels=c("Non-malapportioned", "Malapportioned")) +
  theme_classic() + theme(legend.position="bottom") +
  labs(x = "Year", y = "State Taxes per White Capita")
figA6
#ggsave("figA6.png", width = 5, height = 5, dpi = 300)

# ----------------------------------------------------------------------------
#   Table A4: Taxation Outcomes: Pre-Treatment Parallel Trends, 1835-1844
# ----------------------------------------------------------------------------

# Create trend variable
data$trend <- data$year - 1834
# Subset pre-treatment data
datatrend <- data[which(data$year < 1845),]

# Differences in slope before commodities shock
inc.trend <- plm(st_taxes_per_income_5y ~  malapportion*trend, data = datatrend, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.inc.trend <- coeftest(inc.trend, vcovHC(inc.trend, type = 'HC0', cluster = 'group'))[,2]

pwc.trend <- plm(state_taxes_PWC_5y ~ malapportion*trend, data = datatrend, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.pwc.trend <- coeftest(pwc.trend, vcovHC(pwc.trend, type = 'HC0', cluster = 'group'))[,2]

stargazer(inc.trend, pwc.trend,
          se = list(se.inc.trend, se.pwc.trend), 
          title = "Taxation Outcomes: Pre-Treatment Parallel Trends, 1835-1844",
          dep.var.labels = c("Taxes/Income", "Taxes per White Capita"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes"),
                           c("State FE", "Yes", "Yes"),
                           c("Year FE", "Yes", "Yes"),
                           c("Number of States","12", "12")))


# -----------------------------------------------------------------------------------
#   Table A5: First-Difference Models: Commodity Price Index and Taxation Outcomes
# -----------------------------------------------------------------------------------

fd.inc <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.inc <- coeftest(fd.inc, vcovHC(fd.inc, type = 'HC0', cluster = 'group'))

fd.pwc <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.pwc <- coeftest(fd.pwc, vcovHC(fd.pwc, type = 'HC0', cluster = 'group'))

fd.slv <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.slv <- coeftest(fd.slv, vcovHC(fd.slv, type = 'HC0', cluster = 'group'))

fd.land <- plm(update(base.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.land <- coeftest(fd.land, vcovHC(fd.land, type = 'HC0', cluster = 'group'))

fd.poll <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.poll <- coeftest(fd.poll, vcovHC(fd.poll, type = 'HC0', cluster = 'group'))

stargazer(fd.inc, fd.pwc, fd.slv, fd.land, fd.poll,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index (log) $times$ Malapportionment"),
          title = "First-Difference Models: Commodity Price Index and Taxation Outcomes",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"), 
                           c("Number of States","14", "14", "14", "14", "14")))

# ------------------------------------------------------------------------
#  Table A6: First-Difference Models: Cotton Prices and Taxation Outcomes
# ------------------------------------------------------------------------

fd.inc.c <- plm(update(cotton.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
summary(fd.inc.c)

fd.pwc.c <- plm(update(cotton.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
summary(fd.pwc.c)

fd.slv.c <- plm(update(cotton.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
summary(fd.slv.c)

fd.land.c <- plm(update(cotton.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
summary(fd.land.c)

fd.poll.c <- plm(update(cotton.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
summary(fd.poll.c)

stargazer(fd.inc.c, fd.pwc.c, fd.slv.c, fd.land.c, fd.poll.c,
          keep = "cotton.newOrleans_5y",
          title = "First-Difference Models: Cotton Prices and Taxation Outcomes",
          covariate.labels = c("Cotton Prices (log)", "Cotton Prices $times$ Malapportionment"),
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"), 
                           c("Number of States","14", "14", "14", "14", "14")) )


# --------------------------------------------------
#     Table A7: State Taxes as a Share of Income
# --------------------------------------------------

# specifications
form1 <- formula(. ~ log(index.rel_5y) + malapportion + log(tot_pop+1) + log(urban_pop+1) + log(income_linear))
form2 <- formula(. ~ log(index.rel_5y)*malapportion)

# samples
seceding <- pdata[which(pdata$statefips %!in% c(21,24,29)),]
novotetax <- pdata[which(pdata$statefips %!in% c(37,45)),]
rural <- pdata[which(pdata$statefips %!in% c(24,29,22)),]

# Panel A: Two-Way Fixed Effects Models
# Column 1: w/o interaction
fe.inc.0 <- plm(update(form1, st_taxes_per_income_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.inc.0 <- coeftest(fe.inc.0, vcovHC(fe.inc.0, type = 'HC0', cluster = 'group'))[,2]

# Column 2: w/o controls
fe.inc.1 <- plm(update(form2, st_taxes_per_income_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.inc.1 <- coeftest(fe.inc.1, vcovHC(fe.inc.1, type = 'HC0', cluster = 'group'))[,2]

# Column 3: baseline model
fe.inc.2 <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.inc.2 <- coeftest(fe.inc.2, vcovHC(fe.inc.2, type = 'HC0', cluster = 'group'))[,2]

# Column 4: only seceding states
fe.inc.3 <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = seceding, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.inc.3 <- coeftest(fe.inc.3, vcovHC(fe.inc.3, type = 'HC0', cluster = 'group'))[,2]

# Column 5: excluding states with a vote-tax link
fe.inc.4 <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.inc.4 <- coeftest(fe.inc.4, vcovHC(fe.inc.4, type = 'HC0', cluster = 'group'))[,2]

# Column 6: excluding states with share of urban population
fe.inc.5 <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = rural, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.inc.5 <- coeftest(fe.inc.5, vcovHC(fe.inc.5, type = 'HC0', cluster = 'group'))[,2]

stargazer(fe.inc.0, fe.inc.1, fe.inc.2, fe.inc.3, fe.inc.4, fe.inc.5,
          se = list(se.fe.inc.0, se.fe.inc.1, se.fe.inc.2, se.fe.inc.3, se.fe.inc.4, se.fe.inc.5), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "Two-Way Fixed Effects Models: State Taxes as a Share of Income",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Tax Revenue/Income",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes")) )

# Panel B: First Difference Models
# Column 1: w/o interaction
fd.inc.0 <- plm(update(form1, st_taxes_per_income_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.inc.0)

# Column 2: w/o controls
fd.inc.1 <- plm(update(form2, st_taxes_per_income_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.inc.1)

# Column 3: baseline model
fd.inc.2 <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.inc.2)

# Column 4: only seceding states
fd.inc.3 <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = seceding, index = c("statefips", "year"), model = "fd")
summary(fd.inc.3)

# Column 5: excluding states with a vote-tax link
fd.inc.4 <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "fd") 
summary(fd.inc.4)

# Column 6: excluding states with share of urban population
fd.inc.5 <- plm(update(base.form, st_taxes_per_income_5y ~ .), data = rural, index = c("statefips", "year"), model = "fd") 
summary(fd.inc.5)

stargazer(fd.inc.0, fd.inc.1, fd.inc.2, fd.inc.3, fd.inc.4, fd.inc.5,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "First Difference Models: State Taxes as a Share of Income",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Tax Revenue/Income",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes"),
                           c("States","14", "14", "14", "11", "12", "11") ) )


# ---------------------------------------------------
#    Table A8: State Taxes per White Capita
# ---------------------------------------------------

# Panel A: Two-Way Fixed Effects Models
# Column 1: w/o interaction
fe.pwc.0 <- plm(update(form1, state_taxes_PWC_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.pwc.0 <- coeftest(fe.pwc.0, vcovHC(fe.pwc.0, type = 'HC0', cluster = 'group'))[,2]

# Column 2: w/o controls
fe.pwc.1 <- plm(update(form2, state_taxes_PWC_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.pwc.1 <- coeftest(fe.pwc.1, vcovHC(fe.pwc.1, type = 'HC0', cluster = 'group'))[,2]

# Column 3: baseline model
fe.pwc.2 <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.pwc.2 <- coeftest(fe.pwc.2, vcovHC(fe.pwc.2, type = 'HC0', cluster = 'group'))[,2]

# Column 4: only seceding states
fe.pwc.3 <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = seceding, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.pwc.3 <- coeftest(fe.pwc.3, vcovHC(fe.pwc.3, type = 'HC0', cluster = 'group'))[,2]

# Column 5: excluding states with a vote-tax link
fe.pwc.4 <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.pwc.4 <- coeftest(fe.pwc.4, vcovHC(fe.pwc.4, type = 'HC0', cluster = 'group'))[,2]

# Column 6: excluding states with share of urban population
fe.pwc.5 <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = rural, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.pwc.5 <- coeftest(fe.pwc.5, vcovHC(fe.pwc.5, type = 'HC0', cluster = 'group'))[,2]

stargazer(fe.pwc.0, fe.pwc.1, fe.pwc.2, fe.pwc.3, fe.pwc.4, fe.pwc.5,
          se = list(se.fe.pwc.0, se.fe.pwc.1, se.fe.pwc.2, se.fe.pwc.3, se.fe.pwc.4, se.fe.pwc.5), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "Two-Way Fixed Effects Models: State Taxes per White Capita",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Tax Revenue/White Population",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes")) )

# Panel B: First Difference Models
# Column 1: w/o interaction
fd.pwc.0 <- plm(update(form1, state_taxes_PWC_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.pwc.0)

# Column 2: w/o controls
fd.pwc.1 <- plm(update(form2, state_taxes_PWC_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.pwc.1)

# Column 3: baseline model
fd.pwc.2 <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.pwc.2)

# Column 4: only seceding states
fd.pwc.3 <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = seceding, index = c("statefips", "year"), model = "fd")
summary(fd.pwc.3)

# Column 5: excluding states with a vote-tax link
fd.pwc.4 <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "fd") 
summary(fd.pwc.4)

# Column 6: excluding states with share of urban population
fd.pwc.5 <- plm(update(base.form, state_taxes_PWC_5y ~ .), data = rural, index = c("statefips", "year"), model = "fd") 
summary(fd.pwc.5)

stargazer(fd.pwc.0, fd.pwc.1, fd.pwc.2, fd.pwc.3, fd.pwc.4, fd.pwc.5,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "First Difference Models: State Taxes per White Capita",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Tax Revenue/White Population",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes"),
                           c("States","14", "14", "14", "11", "12", "11") ) )


# ------------------------------------------
#     Table A9: Tax Rate on Slaves (log)
# ------------------------------------------
base.advalorem <- formula(. ~ log(index.rel_5y)*malapportion + log(tot_pop+1) + log(urban_pop+1) + log(income_linear) + advalorem)

# Panel A: Two-Way Fixed Effects Models
# Column 1: w/o interaction
fe.str.0 <- plm(update(form1, log(slaves_tax_rate_5y) ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.str.0 <- coeftest(fe.str.0, vcovHC(fe.str.0, type = 'HC0', cluster = 'group'))[,2]

# Column 2: w/o controls
fe.str.1 <- plm(update(form2, log(slaves_tax_rate_5y) ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.str.1 <- coeftest(fe.str.1, vcovHC(fe.str.1, type = 'HC0', cluster = 'group'))[,2]

# Column 3: baseline model
fe.str.2 <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.str.2 <- coeftest(fe.str.2, vcovHC(fe.str.2, type = 'HC0', cluster = 'group'))[,2]

# Column 4: only seceding states
fe.str.3 <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = seceding, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.str.3 <- coeftest(fe.str.3, vcovHC(fe.str.3, type = 'HC0', cluster = 'group'))[,2]

# Column 5: excluding states with a vote-tax link
fe.str.4 <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = novotetax, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.str.4 <- coeftest(fe.str.4, vcovHC(fe.str.4, type = 'HC0', cluster = 'group'))[,2]

# Column 6: excluding states with share of urban population
fe.str.5 <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = rural, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.str.5 <- coeftest(fe.str.5, vcovHC(fe.str.5, type = 'HC0', cluster = 'group'))[,2]

# Column 7: including dummy for states that had Ad Valorem rather than capitation tax
fe.str.6 <- plm(update(base.advalorem, log(slaves_tax_rate_5y) ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.str.6 <- coeftest(fe.str.6, vcovHC(fe.str.6, type = 'HC0', cluster = 'group'))[,2]


stargazer(fe.str.0, fe.str.1, fe.str.2, fe.str.3, fe.str.4, fe.str.5, fe.str.6,
          se = list(se.fe.str.0, se.fe.str.1, se.fe.str.2, se.fe.str.3, se.fe.str.4, se.fe.str.5, se.fe.str.6), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "Two-Way Fixed Effects Models: Tax Rate on Slaves (log)",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Tax Rate on Slaves (log)",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes")) )

# Panel B: First Difference Models
# Column 1: w/o interaction
fd.str.0 <- plm(update(form1, log(slaves_tax_rate_5y) ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.str.0)

# Column 2: w/o controls
fd.str.1 <- plm(update(form2, log(slaves_tax_rate_5y) ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.str.1)

# Column 3: baseline model
fd.str.2 <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.str.2)

# Column 4: only seceding states
fd.str.3 <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = seceding, index = c("statefips", "year"), model = "fd")
summary(fd.str.3)

# Column 5: excluding states with a vote-tax link
fd.str.4 <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = novotetax, index = c("statefips", "year"), model = "fd") 
summary(fd.str.4)

# Column 6: excluding states with share of urban population
fd.str.5 <- plm(update(base.form, log(slaves_tax_rate_5y) ~ .), data = rural, index = c("statefips", "year"), model = "fd") 
summary(fd.str.5)

# Column 7: including dummy for states that had Ad Valorem rather than capitation tax
fd.str.6 <- plm(update(base.advalorem, log(slaves_tax_rate_5y) ~ .), data = pdata, index = c("statefips", "year"), model = "fd")
se.fd.str.6 <- coeftest(fd.str.6, vcovHC(fd.str.6, type = 'HC0', cluster = 'group'))[,2]


stargazer(fd.str.0, fd.str.1, fd.str.2, fd.str.3, fd.str.4, fd.str.5, fd.str.6,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "First Difference Models: Tax Rate on Slaves (log)",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Tax Rate on Slaves (log)",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes"),
                           c("States","14", "14", "14", "11", "12", "11", "14") ) )

# ----------------------------------------------------
#    Table A10: Average Ad Valorem Tax Rate on Land
# ----------------------------------------------------

# Panel A: Two-Way Fixed Effects Models
# Column 1: w/o interaction
fe.land.0 <- plm(update(form1, tax_rate_land_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.land.0 <- coeftest(fe.land.0, vcovHC(fe.land.0, type = 'HC0', cluster = 'group'))[,2]

# Column 2: w/o controls
fe.land.1 <- plm(update(form2, tax_rate_land_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.land.1 <- coeftest(fe.land.1, vcovHC(fe.land.1, type = 'HC0', cluster = 'group'))[,2]

# Column 3: baseline model
fe.land.2 <- plm(update(base.form, tax_rate_land_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.land.2 <- coeftest(fe.land.2, vcovHC(fe.land.2, type = 'HC0', cluster = 'group'))[,2]

# Column 4: only seceding states
fe.land.3 <- plm(update(base.form, tax_rate_land_5y ~ .), data = seceding, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.land.3 <- coeftest(fe.land.3, vcovHC(fe.land.3, type = 'HC0', cluster = 'group'))[,2]

# Column 5: excluding states with a vote-tax link
fe.land.4 <- plm(update(base.form, tax_rate_land_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.land.4 <- coeftest(fe.land.4, vcovHC(fe.land.4, type = 'HC0', cluster = 'group'))[,2]

# Column 6: excluding states with share of urban population
fe.land.5 <- plm(update(base.form, tax_rate_land_5y ~ .), data = rural, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.land.5 <- coeftest(fe.land.5, vcovHC(fe.land.5, type = 'HC0', cluster = 'group'))[,2]

stargazer(fe.land.0, fe.land.1, fe.land.2, fe.land.3, fe.land.4, fe.land.5, 
          se = list(se.fe.land.0, se.fe.land.1, se.fe.land.2, se.fe.land.3, se.fe.land.4, se.fe.land.5), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "Two-Way Fixed Effects Models: Ad Valorem Tax Rate on Land",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Ad Valorem Tax Rate on Land",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes")) )


# Panel B: First Difference Models
# Column 1: w/o interaction
fd.land.0 <- plm(update(form1, tax_rate_land_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.land.0)

# Column 2: w/o controls
fd.land.1 <- plm(update(form2, tax_rate_land_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.land.1)

# Column 3: baseline model
fd.land.2 <- plm(update(base.form, tax_rate_land_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.land.2)

# Column 4: only seceding states
fd.land.3 <- plm(update(base.form, tax_rate_land_5y ~ .), data = seceding, index = c("statefips", "year"), model = "fd")
summary(fd.land.3)

# Column 5: excluding states with a vote-tax link
fd.land.4 <- plm(update(base.form, tax_rate_land_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "fd") 
summary(fd.land.4)

# Column 6: excluding states with share of urban population
fd.land.5 <- plm(update(base.form, tax_rate_land_5y ~ .), data = rural, index = c("statefips", "year"), model = "fd") 
summary(fd.land.5)


stargazer(fd.land.0, fd.land.1, fd.land.2, fd.land.3, fd.land.4, fd.land.5,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "First Difference Models: Ad Valorem Tax Rate on Land",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Ad Valorem Tax Rate on Land",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes"),
                           c("States","14", "14", "14", "11", "12", "11") ) )

# ----------------------------------------
#    Table A11: Poll Tax per White Male
# ----------------------------------------

# Panel A: Two-Way Fixed Effects Models
# Column 1: w/o interaction
fe.poll.0 <- plm(update(form1, poll_tax_per_wh_male_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.poll.0 <- coeftest(fe.poll.0, vcovHC(fe.poll.0, type = 'HC0', cluster = 'group'))[,2]

# Column 2: w/o controls
fe.poll.1 <- plm(update(form2, poll_tax_per_wh_male_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.poll.1 <- coeftest(fe.poll.1, vcovHC(fe.poll.1, type = 'HC0', cluster = 'group'))[,2]

# Column 3: baseline model
fe.poll.2 <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.poll.2 <- coeftest(fe.poll.2, vcovHC(fe.poll.2, type = 'HC0', cluster = 'group'))[,2]

# Column 4: only seceding states
fe.poll.3 <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = seceding, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.poll.3 <- coeftest(fe.poll.3, vcovHC(fe.poll.3, type = 'HC0', cluster = 'group'))[,2]

# Column 5: excluding states with a vote-tax link
fe.poll.4 <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.poll.4 <- coeftest(fe.poll.4, vcovHC(fe.poll.4, type = 'HC0', cluster = 'group'))[,2]

# Column 6: excluding states with share of urban population
fe.poll.5 <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = rural, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.poll.5 <- coeftest(fe.poll.5, vcovHC(fe.poll.5, type = 'HC0', cluster = 'group'))[,2]

stargazer(fe.poll.0, fe.poll.1, fe.poll.2, fe.poll.3, fe.poll.4, fe.poll.5, 
          se = list(se.fe.poll.0, se.fe.poll.1, se.fe.poll.2, se.fe.poll.3, se.fe.poll.4, se.fe.poll.5), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "Two-Way Fixed Effects Models: Poll Tax per White Male",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Poll Tax Rate",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes")) )

# Panel B: First Difference Models
# Column 1: w/o interaction
fd.poll.0 <- plm(update(form1, poll_tax_per_wh_male_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.poll.0)

# Column 2: w/o controls
fd.poll.1 <- plm(update(form2, poll_tax_per_wh_male_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.poll.1)

# Column 3: baseline model
fd.poll.2 <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.poll.2)

# Column 4: only seceding states
fd.poll.3 <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = seceding, index = c("statefips", "year"), model = "fd")
summary(fd.poll.3)

# Column 5: excluding states with a vote-tax link
fd.poll.4 <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "fd") 
summary(fd.poll.4)

# Column 6: excluding states with share of urban population
fd.poll.5 <- plm(update(base.form, poll_tax_per_wh_male_5y ~ .), data = rural, index = c("statefips", "year"), model = "fd") 
summary(fd.poll.5)


stargazer(fd.poll.0, fd.poll.1, fd.poll.2, fd.poll.3, fd.poll.4, fd.poll.5,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "First Difference Models: Poll Tax per White Male",
          dep.var.labels.include = FALSE,
          dep.var.caption  = "Poll Tax Rate",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes"),
                           c("States","14", "14", "14", "11", "12", "11") ) )


# ----------------------------------------------------------
#  Table A12: Commodity Price Index and Taxation Outcomes with Size of Enslaved Population (log) as Control
# ----------------------------------------------------------
base.form.slv <- formula(. ~ log(index.rel_5y)*malapportion + log(slave_pop+1) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear))

# TWFE
fe.inc.slv <- plm(update(base.form.slv, st_taxes_per_income_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoway")
se.fe.inc.slv <- coeftest(fe.inc.slv, vcovHC(fe.inc.slv, type = 'HC0', cluster = 'group'))[,2]

fe.pwc.slv <- plm(update(base.form.slv, state_taxes_PWC_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoway")
se.fe.pwc.slv <- coeftest(fe.pwc.slv, vcovHC(fe.pwc.slv, type = 'HC0', cluster = 'group'))[,2]

fe.land.slv <- plm(update(base.form.slv, tax_rate_land_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoway")
se.fe.land.slv <- coeftest(fe.land.slv, vcovHC(fe.land.slv, type = 'HC0', cluster = 'group'))[,2]

fe.slv.slv <- plm(update(base.form.slv, log(slaves_tax_rate_5y) ~ .), data = data[which(data$year > 1850),], index = c("statefips", "year"), model = "within", effect = "twoway")
se.fe.slv.slv <- coeftest(fe.slv.slv, vcovHC(fe.slv.slv, type = 'HC0', cluster = 'group'))[,2]

fe.poll.slv <- plm(update(base.form.slv, poll_tax_per_wh_male_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoway")
se.fe.poll.slv <- coeftest(fe.poll.slv, vcovHC(fe.poll.slv, type = 'HC0', cluster = 'group'))[,2]

stargazer(fe.inc.slv, fe.pwc.slv, fe.slv.slv, fe.land.slv, fe.poll.slv,
          keep = c("index.rel_5y", "slave_pop"),
          se = list(se.fe.inc.slv, se.fe.pwc.slv, se.fe.slv.slv, se.fe.land.slv, se.fe.poll.slv), 
          covariate.labels = c("Commodity Index (log)", "Enslaved Population (log)", "Commodity Index $times$ Malapportionment"),
          title = "Two-Way Fixed Effects: Commodity Price Index and Taxation Outcomes with Size of Enslaved Population (log) as Control",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Tax per White Capita", "Tax Rate on Slaves", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("State FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes")))

# FD
fd.inc.slv <- plm(update(base.form.slv, st_taxes_per_income_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd")
se.fd.inc.slv <- coeftest(fd.inc.slv, vcovHC(fd.inc.slv, type = 'HC0', cluster = 'group'))[,2]

fd.pwc.slv <- plm(update(base.form.slv, state_taxes_PWC_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd")
se.fd.pwc.slv <- coeftest(fd.pwc.slv, vcovHC(fd.pwc.slv, type = 'HC0', cluster = 'group'))[,2]

fd.land.slv <- plm(update(base.form.slv, tax_rate_land_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd")
se.fd.land.slv <- coeftest(fd.land.slv, vcovHC(fd.land.slv, type = 'HC0', cluster = 'group'))[,2]

fd.slv.slv <- plm(update(base.form.slv, log(slaves_tax_rate_5y) ~ .), data = data[which(data$year > 1850),], index = c("statefips", "year"), model = "fd")
se.fd.slv.slv <- coeftest(fd.slv.slv, vcovHC(fd.slv.slv, type = 'HC0', cluster = 'group'))[,2]

fd.poll.slv <- plm(update(base.form.slv, poll_tax_per_wh_male_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd")
se.fd.poll.slv <- coeftest(fd.poll.slv, vcovHC(fd.poll.slv, type = 'HC0', cluster = 'group'))[,2]

stargazer(fd.inc.slv, fd.pwc.slv, fd.slv.slv, fd.land.slv, fd.poll.slv,
          keep = c("index.rel_5y", "slave_pop"),
          covariate.labels = c("Commodity Index (log)", "Enslaved Population (log)", "Commodity Index $times$ Malapportionment"),
          title = "First Difference Results: Commodity Price Index and Taxation Outcomes with Size of Enslaved Population (log) as Control",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Tax per White Capita",
                             "Tax Rate on Slaves", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes")))


# ----------------------------------------------------------
#   Table A13: Alternative Moderator: Year of Statehood
# ----------------------------------------------------------
ya.form <- formula(. ~ log(index.rel_5y)*year_admission + log(tot_pop+1) + log(urban_pop+1) + log(income_linear))

# TWFE
fe.ya.inc <- plm(update(ya.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ya.inc <- coeftest(fe.ya.inc, vcovHC(fe.ya.inc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ya.inc <- coeftest(fe.ya.inc, vcovHC(fe.ya.inc, type = 'HC0', cluster = 'group'))

fe.ya.pwc <- plm(update(ya.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ya.pwc <- coeftest(fe.ya.pwc, vcovHC(fe.ya.pwc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ya.pwc <- coeftest(fe.ya.pwc, vcovHC(fe.ya.pwc, type = 'HC0', cluster = 'group'))

fe.ya.slv <- plm(update(ya.form, log(slaves_tax_rate_5y)  ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ya.slv <- coeftest(fe.ya.slv, vcovHC(fe.ya.slv, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ya.slv <- coeftest(fe.ya.slv, vcovHC(fe.ya.slv, type = 'HC0', cluster = 'group'))

fe.ya.land <- plm(update(ya.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ya.land <- coeftest(fe.ya.land, vcovHC(fe.ya.land, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ya.land <- coeftest(fe.ya.land, vcovHC(fe.ya.land, type = 'HC0', cluster = 'group'))

fe.ya.poll <- plm(update(ya.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ya.poll <- coeftest(fe.ya.poll, vcovHC(fe.ya.poll, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ya.poll <- coeftest(fe.ya.poll, vcovHC(fe.ya.poll, type = 'HC0', cluster = 'group'))

stargazer(fe.ya.inc, fe.ya.pwc, fe.ya.slv, fe.ya.land, fe.ya.poll,
          se = list(se.fe.ya.inc, se.fe.ya.pwc, se.fe.ya.slv, se.fe.ya.land, se.fe.ya.poll), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Year of statehood"),
          title = "Two-Way Fixed Effects",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))

# Wild-cluster bootstrap p-values
fe.ya.inc.wild <- cluster.wild.plm(mod=fe.ya.inc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ya.inc.wild)[5]
fe.ya.pwc.wild <- cluster.wild.plm(mod=fe.ya.pwc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ya.pwc.wild)[5]
fe.ya.slv.wild <- cluster.wild.plm(mod=fe.ya.slv, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ya.slv.wild)[5]
fe.ya.land.wild <- cluster.wild.plm(mod=fe.ya.land, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ya.land.wild)[5]
fe.ya.poll.wild <- cluster.wild.plm(mod=fe.ya.poll, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ya.poll.wild)[5]

# First Difference
fd.ya.inc <- plm(update(ya.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ya.inc <- coeftest(fd.ya.inc, vcovHC(fd.ya.inc, type = 'HC0', cluster = 'group'))

fd.ya.pwc <- plm(update(ya.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ya.pwc <- coeftest(fd.ya.pwc, vcovHC(fd.ya.pwc, type = 'HC0', cluster = 'group'))

fd.ya.slv <- plm(update(ya.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ya.slv <- coeftest(fd.ya.slv, vcovHC(fd.ya.slv, type = 'HC0', cluster = 'group'))

fd.ya.land <- plm(update(ya.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ya.land <- coeftest(fd.ya.land, vcovHC(fd.ya.land, type = 'HC0', cluster = 'group'))

fd.ya.poll <- plm(update(ya.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ya.poll <- coeftest(fd.ya.poll, vcovHC(fd.ya.poll, type = 'HC0', cluster = 'group'))

stargazer(fd.ya.inc, fd.ya.pwc, fd.ya.slv, fd.ya.land, fd.ya.poll,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Year of statehood"),
          title = "First Difference",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))

# -------------------------------------------------------------
#   Table A14: Alternative Moderator: Cotton State Indicator
# -------------------------------------------------------------
cot.form <- formula(. ~ log(index.rel_5y)*cotton.state + log(tot_pop+1) + log(urban_pop+1) + log(income_linear))

# TWFE
fe.cot.inc <- plm(update(cot.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.cot.inc <- coeftest(fe.cot.inc, vcovHC(fe.cot.inc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.cot.inc <- coeftest(fe.cot.inc, vcovHC(fe.cot.inc, type = 'HC0', cluster = 'group'))

fe.cot.pwc <- plm(update(cot.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.cot.pwc <- coeftest(fe.cot.pwc, vcovHC(fe.cot.pwc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.cot.pwc <- coeftest(fe.cot.pwc, vcovHC(fe.cot.pwc, type = 'HC0', cluster = 'group'))

fe.cot.slv <- plm(update(cot.form, log(slaves_tax_rate_5y)  ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.cot.slv <- coeftest(fe.cot.slv, vcovHC(fe.cot.slv, type = 'HC0', cluster = 'group'))[,2]
coef.fe.cot.slv <- coeftest(fe.cot.slv, vcovHC(fe.cot.slv, type = 'HC0', cluster = 'group'))

fe.cot.land <- plm(update(cot.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.cot.land <- coeftest(fe.cot.land, vcovHC(fe.cot.land, type = 'HC0', cluster = 'group'))[,2]
coef.fe.cot.land <- coeftest(fe.cot.land, vcovHC(fe.cot.land, type = 'HC0', cluster = 'group'))

fe.cot.poll <- plm(update(cot.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.cot.poll <- coeftest(fe.cot.poll, vcovHC(fe.cot.poll, type = 'HC0', cluster = 'group'))[,2]
coef.fe.cot.poll <- coeftest(fe.cot.poll, vcovHC(fe.cot.poll, type = 'HC0', cluster = 'group'))

stargazer(fe.cot.inc, fe.cot.pwc, fe.cot.slv, fe.cot.land, fe.cot.poll,
          se = list(se.fe.cot.inc, se.fe.cot.pwc, se.fe.cot.slv, se.fe.cot.land, se.fe.cot.poll), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Cotton State"),
          title = "Two-Way Fixed Effects",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))

# Wild-cluster bootstrap p-values
fe.cot.inc.wild <- cluster.wild.plm(mod=fe.cot.inc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.cot.inc.wild)[5]
fe.cot.pwc.wild <- cluster.wild.plm(mod=fe.cot.pwc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.cot.pwc.wild)[5]
fe.cot.slv.wild <- cluster.wild.plm(mod=fe.cot.slv, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.cot.slv.wild)[5]
fe.cot.land.wild <- cluster.wild.plm(mod=fe.cot.land, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.cot.land.wild)[5]
fe.cot.poll.wild <- cluster.wild.plm(mod=fe.cot.poll, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.cot.poll.wild)[5]

# First Difference
fd.cot.inc <- plm(update(cot.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.cot.inc <- coeftest(fd.cot.inc, vcovHC(fd.cot.inc, type = 'HC0', cluster = 'group'))

fd.cot.pwc <- plm(update(cot.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.cot.pwc <- coeftest(fd.cot.pwc, vcovHC(fd.cot.pwc, type = 'HC0', cluster = 'group'))

fd.cot.slv <- plm(update(cot.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.cot.slv <- coeftest(fd.cot.slv, vcovHC(fd.cot.slv, type = 'HC0', cluster = 'group'))

fd.cot.land <- plm(update(cot.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.cot.land <- coeftest(fd.cot.land, vcovHC(fd.cot.land, type = 'HC0', cluster = 'group'))

fd.cot.poll <- plm(update(cot.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.cot.poll <- coeftest(fd.cot.poll, vcovHC(fd.cot.poll, type = 'HC0', cluster = 'group'))

stargazer(fd.cot.inc, fd.cot.pwc, fd.cot.slv, fd.cot.land, fd.cot.poll,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Cotton State"),
          title = "First Difference",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))


# ----------------------------------------------------------------------
#   Table A15: Alternative Moderator: Statehood before 1800 Indicator
# ----------------------------------------------------------------------
s1800.form <- formula(. ~ log(index.rel_5y)*state1800 + log(tot_pop+1) + log(urban_pop+1) + log(income_linear))

# TWFE
fe.1800.inc <- plm(update(s1800.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.1800.inc <- coeftest(fe.1800.inc, vcovHC(fe.1800.inc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.1800.inc <- coeftest(fe.1800.inc, vcovHC(fe.1800.inc, type = 'HC0', cluster = 'group'))

fe.1800.pwc <- plm(update(s1800.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.1800.pwc <- coeftest(fe.1800.pwc, vcovHC(fe.1800.pwc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.1800.pwc <- coeftest(fe.1800.pwc, vcovHC(fe.1800.pwc, type = 'HC0', cluster = 'group'))

fe.1800.slv <- plm(update(s1800.form, log(slaves_tax_rate_5y)  ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.1800.slv <- coeftest(fe.1800.slv, vcovHC(fe.1800.slv, type = 'HC0', cluster = 'group'))[,2]
coef.fe.1800.slv <- coeftest(fe.1800.slv, vcovHC(fe.1800.slv, type = 'HC0', cluster = 'group'))

fe.1800.land <- plm(update(s1800.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.1800.land <- coeftest(fe.1800.land, vcovHC(fe.1800.land, type = 'HC0', cluster = 'group'))[,2]
coef.fe.1800.land <- coeftest(fe.1800.land, vcovHC(fe.1800.land, type = 'HC0', cluster = 'group'))

fe.1800.poll <- plm(update(s1800.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.1800.poll <- coeftest(fe.1800.poll, vcovHC(fe.1800.poll, type = 'HC0', cluster = 'group'))[,2]
coef.fe.1800.poll <- coeftest(fe.1800.poll, vcovHC(fe.1800.poll, type = 'HC0', cluster = 'group'))

stargazer(fe.1800.inc, fe.1800.pwc, fe.1800.slv, fe.1800.land, fe.1800.poll,
          se = list(se.fe.1800.inc, se.fe.1800.pwc, se.fe.1800.slv, se.fe.1800.land, se.fe.1800.poll), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ State before 1800"),
          title = "Two-Way Fixed Effects",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))

# Wild-cluster bootstrap p-values
fe.1800.inc.wild <- cluster.wild.plm(mod=fe.1800.inc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.1800.inc.wild)[5]
fe.1800.pwc.wild <- cluster.wild.plm(mod=fe.1800.pwc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.1800.pwc.wild)[5]
fe.1800.slv.wild <- cluster.wild.plm(mod=fe.1800.slv, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.1800.slv.wild)[5]
fe.1800.land.wild <- cluster.wild.plm(mod=fe.1800.land, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.1800.land.wild)[5]
fe.1800.poll.wild <- cluster.wild.plm(mod=fe.1800.poll, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.1800.poll.wild)[5]

# First Difference
fd.1800.inc <- plm(update(s1800.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.1800.inc <- coeftest(fd.1800.inc, vcovHC(fd.1800.inc, type = 'HC0', cluster = 'group'))

fd.1800.pwc <- plm(update(s1800.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.1800.pwc <- coeftest(fd.1800.pwc, vcovHC(fd.1800.pwc, type = 'HC0', cluster = 'group'))

fd.1800.slv <- plm(update(s1800.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.1800.slv <- coeftest(fd.1800.slv, vcovHC(fd.1800.slv, type = 'HC0', cluster = 'group'))

fd.1800.land <- plm(update(s1800.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.1800.land <- coeftest(fd.1800.land, vcovHC(fd.1800.land, type = 'HC0', cluster = 'group'))

fd.1800.poll <- plm(update(s1800.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.1800.poll <- coeftest(fd.1800.poll, vcovHC(fd.1800.poll, type = 'HC0', cluster = 'group'))

stargazer(fd.1800.inc, fd.1800.pwc, fd.1800.slv, fd.1800.land, fd.1800.poll,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ State before 1800"),
          title = "First Difference",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))

# -----------------------------------------------------------------
#   Table A16: Alternative Moderator: Mississippi River Indicator
# -----------------------------------------------------------------
ms.form <- formula(. ~ log(index.rel_5y)*mississippi + log(tot_pop+1) + log(urban_pop+1) + log(income_linear))

# TWFE
fe.ms.inc <- plm(update(ms.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ms.inc <- coeftest(fe.ms.inc, vcovHC(fe.ms.inc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ms.inc <- coeftest(fe.ms.inc, vcovHC(fe.ms.inc, type = 'HC0', cluster = 'group'))

fe.ms.pwc <- plm(update(ms.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ms.pwc <- coeftest(fe.ms.pwc, vcovHC(fe.ms.pwc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ms.pwc <- coeftest(fe.ms.pwc, vcovHC(fe.ms.pwc, type = 'HC0', cluster = 'group'))

fe.ms.slv <- plm(update(ms.form, log(slaves_tax_rate_5y)  ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ms.slv <- coeftest(fe.ms.slv, vcovHC(fe.ms.slv, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ms.slv <- coeftest(fe.ms.slv, vcovHC(fe.ms.slv, type = 'HC0', cluster = 'group'))

fe.ms.land <- plm(update(ms.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ms.land <- coeftest(fe.ms.land, vcovHC(fe.ms.land, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ms.land <- coeftest(fe.ms.land, vcovHC(fe.ms.land, type = 'HC0', cluster = 'group'))

fe.ms.poll <- plm(update(ms.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.ms.poll <- coeftest(fe.ms.poll, vcovHC(fe.ms.poll, type = 'HC0', cluster = 'group'))[,2]
coef.fe.ms.poll <- coeftest(fe.ms.poll, vcovHC(fe.ms.poll, type = 'HC0', cluster = 'group'))

stargazer(fe.ms.inc, fe.ms.pwc, fe.ms.slv, fe.ms.land, fe.ms.poll,
          se = list(se.fe.ms.inc, se.fe.ms.pwc, se.fe.ms.slv, se.fe.ms.land, se.fe.ms.poll), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ MS River"),
          title = "Two-Way Fixed Effects",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))

# Wild-cluster bootstrap p-values
fe.ms.inc.wild <- cluster.wild.plm(mod=fe.ms.inc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ms.inc.wild)[5]
fe.ms.pwc.wild <- cluster.wild.plm(mod=fe.ms.pwc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ms.pwc.wild)[5]
fe.ms.slv.wild <- cluster.wild.plm(mod=fe.ms.slv, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ms.slv.wild)[5]
fe.ms.land.wild <- cluster.wild.plm(mod=fe.ms.land, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ms.land.wild)[5]
fe.ms.poll.wild <- cluster.wild.plm(mod=fe.ms.poll, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.ms.poll.wild)[5]

# First Difference
fd.ms.inc <- plm(update(ms.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ms.inc <- coeftest(fd.ms.inc, vcovHC(fd.ms.inc, type = 'HC0', cluster = 'group'))

fd.ms.pwc <- plm(update(ms.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ms.pwc <- coeftest(fd.ms.pwc, vcovHC(fd.ms.pwc, type = 'HC0', cluster = 'group'))

fd.ms.slv <- plm(update(ms.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ms.slv <- coeftest(fd.ms.slv, vcovHC(fd.ms.slv, type = 'HC0', cluster = 'group'))

fd.ms.land <- plm(update(ms.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ms.land <- coeftest(fd.ms.land, vcovHC(fd.ms.land, type = 'HC0', cluster = 'group'))

fd.ms.poll <- plm(update(ms.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.ms.poll <- coeftest(fd.ms.poll, vcovHC(fd.ms.poll, type = 'HC0', cluster = 'group'))

stargazer(fd.ms.inc, fd.ms.pwc, fd.ms.slv, fd.ms.land, fd.ms.poll,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ MS River"),
          title = "First Difference",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))

# ------------------------------------------------------------
#  Table A17: Alternative Moderator: Coastal State Indicator
# ------------------------------------------------------------
coast.form <- formula(. ~ log(index.rel_5y)*coastal + log(tot_pop+1) + log(urban_pop+1) + log(income_linear))

# TWFE
fe.coast.inc <- plm(update(coast.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.coast.inc <- coeftest(fe.coast.inc, vcovHC(fe.coast.inc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.coast.inc <- coeftest(fe.coast.inc, vcovHC(fe.coast.inc, type = 'HC0', cluster = 'group'))

fe.coast.pwc <- plm(update(coast.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.coast.pwc <- coeftest(fe.coast.pwc, vcovHC(fe.coast.pwc, type = 'HC0', cluster = 'group'))[,2]
coef.fe.coast.pwc <- coeftest(fe.coast.pwc, vcovHC(fe.coast.pwc, type = 'HC0', cluster = 'group'))

fe.coast.slv <- plm(update(coast.form, log(slaves_tax_rate_5y)  ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.coast.slv <- coeftest(fe.coast.slv, vcovHC(fe.coast.slv, type = 'HC0', cluster = 'group'))[,2]
coef.fe.coast.slv <- coeftest(fe.coast.slv, vcovHC(fe.coast.slv, type = 'HC0', cluster = 'group'))

fe.coast.land <- plm(update(coast.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.coast.land <- coeftest(fe.coast.land, vcovHC(fe.coast.land, type = 'HC0', cluster = 'group'))[,2]
coef.fe.coast.land <- coeftest(fe.coast.land, vcovHC(fe.coast.land, type = 'HC0', cluster = 'group'))

fe.coast.poll <- plm(update(coast.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="within", effect = "twoways")
se.fe.coast.poll <- coeftest(fe.coast.poll, vcovHC(fe.coast.poll, type = 'HC0', cluster = 'group'))[,2]
coef.fe.coast.poll <- coeftest(fe.coast.poll, vcovHC(fe.coast.poll, type = 'HC0', cluster = 'group'))

stargazer(fe.coast.inc, fe.coast.pwc, fe.coast.slv, fe.coast.land, fe.coast.poll,
          se = list(se.fe.coast.inc, se.fe.coast.pwc, se.fe.coast.slv, se.fe.coast.land, se.fe.coast.poll), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Coastal State"),
          title = "Two-Way Fixed Effects",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on Land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))

# Wild-cluster bootstrap p-values
fe.coast.inc.wild <- cluster.wild.plm(mod=fe.coast.inc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.coast.inc.wild)[5]
fe.coast.pwc.wild <- cluster.wild.plm(mod=fe.coast.pwc, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.coast.pwc.wild)[5]
fe.coast.slv.wild <- cluster.wild.plm(mod=fe.coast.slv, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.coast.slv.wild)[5]
fe.coast.land.wild <- cluster.wild.plm(mod=fe.coast.land, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.coast.land.wild)[5]
fe.coast.poll.wild <- cluster.wild.plm(mod=fe.coast.poll, dat=pdata, cluster="group", ci.level = 0.95, boot.reps = 1000, report = TRUE, prog.bar = TRUE)[1]
unlist(fe.coast.poll.wild)[5]

# First Difference
fd.coast.inc <- plm(update(coast.form, st_taxes_per_income_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.coast.inc <- coeftest(fd.coast.inc, vcovHC(fd.coast.inc, type = 'HC0', cluster = 'group'))

fd.coast.pwc <- plm(update(coast.form, state_taxes_PWC_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.coast.pwc <- coeftest(fd.coast.pwc, vcovHC(fd.coast.pwc, type = 'HC0', cluster = 'group'))

fd.coast.slv <- plm(update(coast.form, log(slaves_tax_rate_5y) ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.coast.slv <- coeftest(fd.coast.slv, vcovHC(fd.coast.slv, type = 'HC0', cluster = 'group'))

fd.coast.land <- plm(update(coast.form, tax_rate_land_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.coast.land <- coeftest(fd.coast.land, vcovHC(fd.coast.land, type = 'HC0', cluster = 'group'))

fd.coast.poll <- plm(update(coast.form, poll_tax_per_wh_male_5y ~ .), data = pdata, index=c("stateabbrev", "year"), model="fd")
coef.fd.coast.poll <- coeftest(fd.coast.poll, vcovHC(fd.coast.poll, type = 'HC0', cluster = 'group'))

stargazer(fd.coast.inc, fd.coast.pwc, fd.coast.slv, fd.coast.land, fd.coast.poll,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Coastal State"),
          title = "First Difference",
          dep.var.labels.include = FALSE,
          column.labels = c("Taxes/Income", "Taxes per White Capita", "Tax Rate on Slaves (log)", "Tax Rate on land", "Poll Tax Rate"),
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("States","14", "14", "14", "14", "14")))

# --------------------------------------------------------
#  Figure A7: Coefficient Estimates of Interaction Terms
# --------------------------------------------------------

# ------------------------------------------
# Figure A7 (a) Coefficient Plot TWFE models
cplot.fe <- rbind(coef.fe.ya.inc, coef.fe.ya.pwc, coef.fe.ya.slv, coef.fe.ya.land, coef.fe.ya.poll)
cplot.fe <- as.data.frame(cplot.fe)
cplot.fe$model <- NA
cplot.fe$model[1:5] <- "Taxes/Income"
cplot.fe$model[6:10] <- "Taxes/White pop"
cplot.fe$model[11:15] <- "Tax rate on slaves"
cplot.fe$model[16:20] <- "Tax rate on land"
cplot.fe$model[21:25] <- "Poll tax rate"
colnames(cplot.fe) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot.fe <- cplot.fe[c(5,10,15,20, 25),]
cplot.fe$term <- "Index:Year Admission"

cplot2.fe <- rbind(coef.fe.cot.inc, coef.fe.cot.pwc, coef.fe.cot.slv, coef.fe.cot.land, coef.fe.cot.poll)
cplot2.fe <- as.data.frame(cplot2.fe)
cplot2.fe <- cplot2.fe[c(5,10,15,20, 25),]
cplot2.fe$model <- NA
cplot2.fe$model[1] <- "Taxes/Income"
cplot2.fe$model[2] <- "Taxes/White pop"
cplot2.fe$model[3] <- "Tax rate on slaves"
cplot2.fe$model[4] <- "Tax rate on land"
cplot2.fe$model[5] <- "Poll tax rate"
colnames(cplot2.fe) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot2.fe$term <- "Index:Cotton State"

cplot3.fe <- rbind(coef.fe.1800.inc, coef.fe.1800.pwc, coef.fe.1800.slv, coef.fe.1800.land, coef.fe.1800.poll)
cplot3.fe <- as.data.frame(cplot3.fe)
cplot3.fe <- cplot3.fe[c(5,10,15,20,25),]
cplot3.fe$model <- NA
cplot3.fe$model[1] <- "Taxes/Income"
cplot3.fe$model[2] <- "Taxes/White pop"
cplot3.fe$model[3] <- "Tax rate on slaves"
cplot3.fe$model[4] <- "Tax rate on land"
cplot3.fe$model[5] <- "Poll tax rate"
colnames(cplot3.fe) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot3.fe$term <- "Index:State before 1800"

cplot5.fe <- rbind(coef.fe.ms.inc, coef.fe.ms.pwc, coef.fe.ms.slv, coef.fe.ms.land, coef.fe.ms.poll)
cplot5.fe <- as.data.frame(cplot5.fe)
cplot5.fe <- cplot5.fe[c(5,10,15,20,25),]
cplot5.fe$model <- NA
cplot5.fe$model[1] <- "Taxes/Income"
cplot5.fe$model[2] <- "Taxes/White pop"
cplot5.fe$model[3] <- "Tax rate on slaves"
cplot5.fe$model[4] <- "Tax rate on land"
cplot5.fe$model[5] <- "Poll tax rate"
colnames(cplot5.fe) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot5.fe$term <- "Index:Mississippi River"

cplot4.fe <- rbind(coef.tw.inc, coef.tw.pwc, coef.tw.slv, coef.tw.land, coef.tw.poll)
cplot4.fe <- as.data.frame(cplot4.fe)
cplot4.fe <- cplot4.fe[c(5,10,15,20,25),]
cplot4.fe$model <- NA
cplot4.fe$model[1] <- "Taxes/Income"
cplot4.fe$model[2] <- "Taxes/White pop"
cplot4.fe$model[3] <- "Tax rate on slaves"
cplot4.fe$model[4] <- "Tax rate on land"
cplot4.fe$model[5] <- "Poll tax rate"
colnames(cplot4.fe) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot4.fe$term <- "Index:Malapportion"

cplot6.fe <- rbind(coef.fe.coast.inc, coef.fe.coast.pwc, coef.fe.coast.slv, coef.fe.coast.land, coef.fe.coast.poll)
cplot6.fe <- as.data.frame(cplot6.fe)
cplot6.fe <- cplot6.fe[c(5,10,15,20,25),]
cplot6.fe$model <- NA
cplot6.fe$model[1] <- "Taxes/Income"
cplot6.fe$model[2] <- "Taxes/White pop"
cplot6.fe$model[3] <- "Tax rate on slaves"
cplot6.fe$model[4] <- "Tax rate on land"
cplot6.fe$model[5] <- "Poll tax rate"
colnames(cplot6.fe) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot6.fe$term <- "Index:Coastal"

coef.plot.fe <- rbind(cplot.fe, cplot3.fe, cplot2.fe, cplot5.fe, cplot6.fe, cplot4.fe)
names(coef.plot.fe)[5:6] <- c("term", "model")

library(dotwhisker)
coefplot.fe.index <- small_multiple(coef.plot.fe) +
  theme_bw() + 
  ylab("Coefficient Estimate") +
  xlab("Interaction Term") +
  geom_hline(yintercept = 0, colour = "grey60", linetype = 2) +
  theme(axis.text.x  = element_text(angle = 45, hjust = 1),
        text = element_text(size=15),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none")
coefplot.fe.index
#ggsave("coefplot.fe.index.png", width = 5, height = 10, dpi = 300)

# ----------------------------------------
# Figure A7 (b) Coefficient Plot FD models
cplot <- rbind(coef.fd.ya.inc, coef.fd.ya.pwc, coef.fd.ya.slv, coef.fd.ya.land, coef.fd.ya.poll)
cplot <- as.data.frame(cplot)
cplot$model <- NA
cplot$model[1:6] <- "Taxes/Income"
cplot$model[7:12] <- "Taxes/White pop"
cplot$model[13:18] <- "Tax rate on slaves"
cplot$model[19:24] <- "Tax rate on land"
cplot$model[25:30] <- "Poll tax rate"
colnames(cplot) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot <- cplot[c(6,12,18,24,30),]
cplot$term <- "Index:Year Admission"

cplot2 <- rbind(coef.fd.cot.inc, coef.fd.cot.pwc, coef.fd.cot.slv, coef.fd.cot.land, coef.fd.cot.poll)
cplot2 <- as.data.frame(cplot2)
cplot2 <- cplot2[c(6,12,18,24, 30),]
cplot2$model <- NA
cplot2$model[1] <- "Taxes/Income"
cplot2$model[2] <- "Taxes/White pop"
cplot2$model[3] <- "Tax rate on slaves"
cplot2$model[4] <- "Tax rate on land"
cplot2$model[5] <- "Poll tax rate"
colnames(cplot2) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot2$term <- "Index:Cotton State"

cplot3 <- rbind(coef.fd.1800.inc, coef.fd.1800.pwc, coef.fd.1800.slv, coef.fd.1800.land, coef.fd.1800.poll)
cplot3 <- as.data.frame(cplot3)
cplot3 <- cplot3[c(6,12,18,24, 30),]
cplot3$model <- NA
cplot3$model[1] <- "Taxes/Income"
cplot3$model[2] <- "Taxes/White pop"
cplot3$model[3] <- "Tax rate on slaves"
cplot3$model[4] <- "Tax rate on land"
cplot3$model[5] <- "Poll tax rate"
colnames(cplot3) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot3$term <- "Index:State before 1800"

cplot5 <- rbind(coef.fd.ms.inc, coef.fd.ms.pwc, coef.fd.ms.slv, coef.fd.ms.land, coef.fd.ms.poll)
cplot5 <- as.data.frame(cplot5)
cplot5 <- cplot5[c(6,12,18,24, 30),]
cplot5$model <- NA
cplot5$model[1] <- "Taxes/Income"
cplot5$model[2] <- "Taxes/White pop"
cplot5$model[3] <- "Tax rate on slaves"
cplot5$model[4] <- "Tax rate on land"
cplot5$model[5] <- "Poll tax rate"
colnames(cplot5) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot5$term <- "Index:Mississippi River"

cplot4 <- rbind(coef.fd.inc, coef.fd.pwc, coef.fd.slv, coef.fd.land, coef.fd.poll)
cplot4 <- as.data.frame(cplot4)
cplot4 <- cplot4[c(6,12,18,24, 30),]
cplot4$model <- NA
cplot4$model[1] <- "Taxes/Income"
cplot4$model[2] <- "Taxes/White pop"
cplot4$model[3] <- "Tax rate on slaves"
cplot4$model[4] <- "Tax rate on land"
cplot4$model[5] <- "Poll tax rate"
colnames(cplot4) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot4$term <- "Index:Malapportion"

cplot6 <- rbind(coef.fd.coast.inc, coef.fd.coast.pwc, coef.fd.coast.slv, coef.fd.coast.land, coef.fd.coast.poll)
cplot6 <- as.data.frame(cplot6)
cplot6 <- cplot6[c(6,12,18,24,30),]
cplot6$model <- NA
cplot6$model[1] <- "Taxes/Income"
cplot6$model[2] <- "Taxes/White pop"
cplot6$model[3] <- "Tax rate on slaves"
cplot6$model[4] <- "Tax rate on land"
cplot6$model[5] <- "Poll tax rate"
colnames(cplot6) <- c("estimate", "std.error", "statistic", "p.value", "model")
cplot6$term <- "Index:Coastal"


coef.plot <- rbind(cplot, cplot3, cplot2, cplot5, cplot6, cplot4)
names(coef.plot)[5:6] <- c("term", "model")

coefplot.fd.index <- small_multiple(coef.plot) +
  theme_bw() + 
  ylab("Coefficient Estimate") +
  xlab("Interaction Term") +
  geom_hline(yintercept = 0, colour = "grey60", linetype = 2) +
  theme(axis.text.x  = element_text(angle = 45, hjust = 1),
        text = element_text(size=15),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "none") +
  scale_y_continuous(breaks = c(-.5,0,.5,1))
coefplot.fd.index
#ggsave("coefplot.fd.index.png", width = 5, height = 10, dpi = 300)


# ---------------------------------------------------------------
#  Table A18: Total Poll Tax Revenues per White Male Capita ($)
# ---------------------------------------------------------------

# Two-Way Fixed Effects Models
# Column 1: w/o interaction
fe.pollrev.0 <- plm(update(form1, tot_poll_tax_rev_wh_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.pollrev.0 <- coeftest(fe.pollrev.0, vcovHC(fe.pollrev.0, type = 'HC0', cluster = 'group'))[,2]

# Column 2: w/o controls
fe.pollrev.1 <- plm(update(form2, tot_poll_tax_rev_wh_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.pollrev.1 <- coeftest(fe.pollrev.1, vcovHC(fe.pollrev.1, type = 'HC0', cluster = 'group'))[,2]

# Column 3: baseline model
fe.pollrev.2 <- plm(update(base.form, tot_poll_tax_rev_wh_5y ~ .), data = pdata, index = c("statefips", "year"), model = "within", effect = "twoways") 
se.fe.pollrev.2 <- coeftest(fe.pollrev.2, vcovHC(fe.pollrev.2, type = 'HC0', cluster = 'group'))[,2]

# Column 4: only seceding states
fe.pollrev.3 <- plm(update(base.form, tot_poll_tax_rev_wh_5y ~ .), data = seceding, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.pollrev.3 <- coeftest(fe.pollrev.3, vcovHC(fe.pollrev.3, type = 'HC0', cluster = 'group'))[,2]

# Column 5: excluding states with a vote-tax link
fe.pollrev.4 <- plm(update(base.form, tot_poll_tax_rev_wh_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.pollrev.4 <- coeftest(fe.pollrev.4, vcovHC(fe.pollrev.4, type = 'HC0', cluster = 'group'))[,2]

# Column 6: excluding states with share of urban population
fe.pollrev.5 <- plm(update(base.form, tot_poll_tax_rev_wh_5y ~ .), data = rural, index = c("statefips", "year"), model = "within", effect = "twoways")
se.fe.pollrev.5 <- coeftest(fe.pollrev.5, vcovHC(fe.pollrev.5, type = 'HC0', cluster = 'group'))[,2]

stargazer(fe.pollrev.0, fe.pollrev.1, fe.pollrev.2, fe.pollrev.3, fe.pollrev.4, fe.pollrev.5, 
          se = list(se.fe.pollrev.0, se.fe.pollrev.1, se.fe.pollrev.2, se.fe.pollrev.3, se.fe.pollrev.4, se.fe.pollrev.5), 
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "Two-Way Fixed Effects Models: Total Poll Tax Revenues per White Male Capita ($)",
          dep.var.labels.include = FALSE,
          dep.var.caption = "Poll Tax Revenue/White Male Population",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes")) )

# First-Difference Models
# Column 1: w/o interaction
fd.pollrev.0 <- plm(update(form1, tot_poll_tax_rev_wh_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.pollrev.0)

# Column 2: w/o controls
fd.pollrev.1 <- plm(update(form2, tot_poll_tax_rev_wh_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.pollrev.1)

# Column 3: baseline model
fd.pollrev.2 <- plm(update(base.form, tot_poll_tax_rev_wh_5y ~ .), data = pdata, index = c("statefips", "year"), model = "fd") 
summary(fd.pollrev.2)

# Column 4: only seceding states
fd.pollrev.3 <- plm(update(base.form, tot_poll_tax_rev_wh_5y ~ .), data = seceding, index = c("statefips", "year"), model = "fd")
summary(fd.pollrev.3)

# Column 5: excluding states with a vote-tax link
fd.pollrev.4 <- plm(update(base.form, tot_poll_tax_rev_wh_5y ~ .), data = novotetax, index = c("statefips", "year"), model = "fd") 
summary(fd.pollrev.4)

# Column 6: excluding states with share of urban population
fd.pollrev.5 <- plm(update(base.form, tot_poll_tax_rev_wh_5y ~ .), data = rural, index = c("statefips", "year"), model = "fd") 
summary(fd.pollrev.5)

stargazer(fd.pollrev.0, fd.pollrev.1, fd.pollrev.2, fd.pollrev.3, fd.pollrev.4, fd.pollrev.5,
          keep = "index.rel_5y",
          covariate.labels = c("Commodity Index (log)", "Commodity Index $times$ Malapportionment"),
          title = "First Difference Models: Total Poll Tax Revenues per White Male Capita ($)",
          dep.var.labels.include = FALSE,
          dep.var.caption = "Poll Tax Revenue/White Male Population",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "No", "Yes", "Yes", "Yes", "Yes"), 
                           c("Only seceding states","No", "No", "No", "Yes", "No", "No"),
                           c("States w/o vote-tax link","No", "No", "No", "No", "Yes", "No"),
                           c("Without urbanized states","No", "No", "No", "No", "No", "Yes"),
                           c("States","14", "14", "14", "11", "12", "11") ) )

# -----------------------------------------------------------------------------------
# Table A19: State Taxes PWC and Tax Rates on Slaves with 1-, 2-, 3-, and 4-Year Lags
# -----------------------------------------------------------------------------------

# Lagged independent variables
lags <- pdata %>% 
  group_by(statefips) %>% 
  mutate(lag.rate1 = dplyr::lag(slaves_tax_rate, n = 1, default = NA),
         lag.rate2 = dplyr::lag(slaves_tax_rate, n = 2, default = NA),
         lag.rate3 = dplyr::lag(slaves_tax_rate, n = 3, default = NA),
         lag.rate4 = dplyr::lag(slaves_tax_rate, n = 4, default = NA),
         lag.rate5 = dplyr::lag(slaves_tax_rate, n = 5, default = NA)) %>% 
  mutate(lag.land1 = dplyr::lag(tax_rate_land, n = 1, default = NA),
         lag.land2 = dplyr::lag(tax_rate_land, n = 2, default = NA),
         lag.land3 = dplyr::lag(tax_rate_land, n = 3, default = NA), 
         lag.land4 = dplyr::lag(tax_rate_land, n = 4, default = NA),
         lag.land5 = dplyr::lag(tax_rate_land, n = 5, default = NA))

# TWFE models
fe.slv.rt0 <- plm(state_taxes_PWC ~ log(slaves_tax_rate) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.slv.rt0 <- coeftest(fe.slv.rt0, vcovHC(fe.slv.rt0, type = 'HC0', cluster = 'group'))[,2]

fe.slv.rt1 <- plm(state_taxes_PWC ~ log(lag.rate1) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.slv.rt1 <- coeftest(fe.slv.rt1, vcovHC(fe.slv.rt1, type = 'HC0', cluster = 'group'))[,2]

fe.slv.rt2 <- plm(state_taxes_PWC ~ log(lag.rate2) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.slv.rt2 <- coeftest(fe.slv.rt2, vcovHC(fe.slv.rt2, type = 'HC0', cluster = 'group'))[,2]

fe.slv.rt3 <- plm(state_taxes_PWC ~ log(lag.rate3) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.slv.rt3 <- coeftest(fe.slv.rt3, vcovHC(fe.slv.rt3, type = 'HC0', cluster = 'group'))[,2]

fe.slv.rt4 <- plm(state_taxes_PWC ~ log(lag.rate4) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.slv.rt4 <- coeftest(fe.slv.rt4, vcovHC(fe.slv.rt4, type = 'HC0', cluster = 'group'))[,2]

stargazer(fe.slv.rt0, fe.slv.rt1, fe.slv.rt2, fe.slv.rt3, fe.slv.rt4,
          se = list(se.fe.slv.rt0, se.fe.slv.rt1, se.fe.slv.rt2, se.fe.slv.rt3, se.fe.slv.rt4), 
          keep = c("slaves_tax_rate", "lag.rate1", "lag.rate2", "lag.rate3", "lag.rate4"),
          covariate.labels = c("Tax rate on slaves", "Tax rate on slaves (1-year lag)","Tax rate on slaves (2-year lag)", "Tax rate on slaves (3-year lag)", "Tax rate on slaves (4-year lag)"),
          title = "State Taxes per White Capita and Tax Rates on Slaves with 1-,2-,3-, and 4-Year Lags",
          dep.var.labels = "Tax Revenue/White Population",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes")))

# FD models
fd.slv.rt0 <- plm(state_taxes_PWC ~ log(slaves_tax_rate) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

fd.slv.rt1 <- plm(state_taxes_PWC ~ log(lag.rate1) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

fd.slv.rt2 <- plm(state_taxes_PWC ~ log(lag.rate2) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

fd.slv.rt3 <- plm(state_taxes_PWC ~ log(lag.rate3) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

fd.slv.rt4 <- plm(state_taxes_PWC ~ log(lag.rate4) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

stargazer(fd.slv.rt0, fd.slv.rt1, fd.slv.rt2, fd.slv.rt3, fd.slv.rt4,
          title = "State Taxes per White Capita and Tax Rates on Slaves with 1-,2-,3-, and 4-Year Lags",
          keep = c("slaves_tax_rate", "lag.rate1", "lag.rate2", "lag.rate3", "lag.rate4"),
          covariate.labels = c("Tax rate on slaves", "Tax rate on slaves (1-year lag)","Tax rate on slaves (2-year lag)", "Tax rate on slaves (3-year lag)", "Tax rate on slaves (4-year lag)"),
          dep.var.labels = "Tax Revenue/White Population",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes")))

# --------------------------------------------------------------------------------
#  Table A20: State Taxes PWC and Land Tax Rate with 1-, 2-, 3-, and 4-Year Lags
# --------------------------------------------------------------------------------

# TWFE models
fe.land.rt0 <- plm(state_taxes_PWC ~ tax_rate_land + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.land.rt0 <- coeftest(fe.land.rt0, vcovHC(fe.land.rt0, type = 'HC0', cluster = 'group'))[,2]

fe.land.rt1 <- plm(state_taxes_PWC ~ log(lag.land1) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.land.rt1 <- coeftest(fe.land.rt1, vcovHC(fe.land.rt1, type = 'HC0', cluster = 'group'))[,2]

fe.land.rt2 <- plm(state_taxes_PWC ~ log(lag.land2) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.land.rt2 <- coeftest(fe.land.rt2, vcovHC(fe.land.rt2, type = 'HC0', cluster = 'group'))[,2]

fe.land.rt3 <- plm(state_taxes_PWC ~ log(lag.land3) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.land.rt3 <- coeftest(fe.land.rt3, vcovHC(fe.land.rt3, type = 'HC0', cluster = 'group'))[,2]

fe.land.rt4 <- plm(state_taxes_PWC ~ log(lag.land4) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index=c("stateabbrev", "year"), model="within", effect = "twoways") 
se.fe.land.rt4 <- coeftest(fe.land.rt4, vcovHC(fe.land.rt4, type = 'HC0', cluster = 'group'))[,2]

stargazer(fe.land.rt0, fe.land.rt1, fe.land.rt2, fe.land.rt3, fe.land.rt4,
          keep = c("tax_rate_land", "lag.land1", "lag.land2", "lag.land3", "lag.land4"),
          covariate.labels = c("Tax rate on land", "Tax rate on land (1-year lag)","Tax rate on land (2-year lag)", "Tax rate on land (3-year lag)", "Tax rate on land (4-year lag)"),
          se = list(se.fe.land.rt0, se.fe.land.rt1, se.fe.land.rt2, se.fe.land.rt3, se.fe.land.rt4), 
          title = "State Taxes per White Capita and Ad Valorem Land Tax with 1-,2-,3-, and 4-Year Lags",
          dep.var.labels = "Tax Revenue/White Population",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("State FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE","Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes")))

# FD models
fd.land.rt0 <- plm(state_taxes_PWC ~ tax_rate_land + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

fd.land.rt1 <- plm(state_taxes_PWC ~ log(lag.land1) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

fd.land.rt2 <- plm(state_taxes_PWC ~ log(lag.land2) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

fd.land.rt3 <- plm(state_taxes_PWC ~ log(lag.land3) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

fd.land.rt4 <- plm(state_taxes_PWC ~ log(lag.land4) + log(tot_pop+1) + log(urban_pop+1) + log(income_linear), data = lags, index = c("statefips", "year"), model = "fd") 

stargazer(fd.land.rt0, fd.land.rt1, fd.land.rt2, fd.land.rt3, fd.land.rt4,
          keep = c("tax_rate_land", "lag.land1", "lag.land2", "lag.land3", "lag.land4"),
          covariate.labels = c("Tax rate on land", "Tax rate on land (1-year lag)","Tax rate on land (2-year lag)", "Tax rate on land (3-year lag)", "Tax rate on land (4-year lag)"),
          title = "State Taxes per White Capita and Ad Valorem Land Tax with 1-,2-,3-, 4- and 5-Year Lags",
          dep.var.labels = "Tax Revenue/White Population",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("Additional covariates","Yes", "Yes", "Yes", "Yes", "Yes")))

# --------------------------------------------------------------------------------
# Table A22: Railroad Miles and Malapportionment, Urbanization, Manufacturing, and Enslaved Share of the Population
# --------------------------------------------------------------------------------

# Railroad Miles/Income
rr.inc1 <- lm(log(rr_income+1) ~ malapportion + slave.pc, data = data)

rr.inc2 <- lm(log(rr_income+1) ~ malapportion + urban.pc, data = data)

rr.inc3 <- lm(log(rr_income+1) ~ malapportion + log(manuf.value.product), data = data)

rr.inc4 <- lm(log(rr_income+1) ~ malapportion + slave.pc + urban.pc + log(manuf.value.product), data = data)

stargazer(rr.inc1, rr.inc2, rr.inc3, rr.inc4,
          covariate.labels = c("Malapportionment", "Enslaved Population (%)", "Urban Population (%)", "Value of Manufacturing Products (log)"),
          title = "Railroad Miles and Malapportionment, Urbanization, Manufacturing, and Enslaved Share of the Population",
          dep.var.labels = "Log of Railroad Miles/Income (Hundred $)",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"))

# Railroad Miles per White Capita (log)
rr.pwc1 <- lm(log(rr_pwc+1) ~ malapportion + slave.pc, data = data)

rr.pwc2 <- lm(log(rr_pwc+1) ~ malapportion + urban.pc, data = data)

rr.pwc3 <- lm(log(rr_pwc+1) ~ malapportion + log(manuf.value.product), data = data)

rr.pwc4 <- lm(log(rr_pwc+1) ~ malapportion + slave.pc + urban.pc + log(manuf.value.product), data = data)

stargazer(rr.pwc1, rr.pwc2, rr.pwc3, rr.pwc4,
          title = "Railroad Miles and Malapportionment, Urbanization, Manufacturing, and Enslaved Share of the Population",
          covariate.labels = c("Malapportionment", "Enslaved Population (%)", "Urban Population (%)", "Value of Manufacturing Products (log)"),
          dep.var.labels = "Log of Railroad Miles per White Capita (Hundred)",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"))

# ------------------------------------------------------------------------------
#     Table A23: County Enslaved Population Share and Miles of Railroads
# ------------------------------------------------------------------------------
rrc <- read.csv("rr_counties.csv")
rrc$slave_sha_1860 <- rrc$slaves_1860/rrc$pop_1860

a1 <- lm(ln_railroads_1860 ~ slave_sha_1860 + pop_1860 + area + ln_sbnr_1860 + log_distance_atlantic + log_distance_gulf + as.factor(stateabbrev), data = rrc[which(rrc$malapportion == 1),])
cl.se.a1 <- coeftest(a1, vcov=function(x) cluster.vcov(x, rrc[which(rrc$malapportion == 1), "stateabbrev"]))[,2]

a0 <- lm(ln_railroads_1860 ~ slave_sha_1860 + pop_1860 + area + ln_sbnr_1860 + log_distance_atlantic + log_distance_gulf + as.factor(stateabbrev), data = rrc[which(rrc$malapportion == 0),])
cl.se.a0 <- coeftest(a0, vcov=function(x) cluster.vcov(x, rrc[which(rrc$malapportion == 0), "stateabbrev"]))[,2]

stargazer(a1, a0,
          se = list(cl.se.a1, cl.se.a0), 
          keep= c("slave_sha_1860"),
          covariate.labels = c("Enslaved Population (%)"),
          title = "County Enslaved-Population Share and Railroad Miles, 1860",
          dep.var.labels = "County Railroad Miles",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("States", "7", "7"),
                           c("State FE", "Yes", "Yes"),
                           c("Additional covariates", "Yes", "Yes")) )

# -----------------------------------------------------
#  Figure B1: Cotton Prices and Commodity Price Index
# -----------------------------------------------------

# Figure B1 (a): Cotton Prices ($)
figB1a <- ggplot(data, aes(x = year)) +
  geom_line(aes(y = cotton.newOrleans_5y)) +
  labs(y = "Price in nominal dollars ($)", x = "Year", colour = "") + 
  scale_colour_manual(values = "black") +
  theme_bw() + xlim(c(1840, 1860)) + ylim(c(6, 12)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
figB1a
ggsave("figB1a.png", width = 5, height = 5, dpi = 300)

# Figure B1 (b):  Commodity Price Index
figB1b_data <- data %>% 
  group_by(year) %>% 
  summarise(mean_group = mean(index.rel_5y, na.rm = TRUE))
figB1b <- ggplot(figB1b_data, aes(x = year, y = mean_group)) +
  geom_line(aes(y = mean_group)) +
  scale_colour_manual(values = "black") + 
  labs(x = "Year", y = "Average Commodity Price Index") +
  theme_bw() + xlim(c(1840, 1860)) + ylim(c(16.5, 31.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
figB1b
ggsave("figB1b.png", width = 5, height = 5, dpi = 300)

# ------------------------------------
#   Figure B2: Average Slave Prices
# ------------------------------------

figB2 <- ggplot(data, aes(x = year)) +
  geom_line(aes(y = slave.price_5y, colour = "")) +
  scale_colour_manual(values = "black") +
  labs(y = "Price in nominal dollars ($)", x = "Year", colour = "") + 
  xlim(c(1840, 1860)) + ylim(c(200, 800)) +
  theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
figB2
ggsave("figB2.png", width = 5, height = 5, dpi = 300)

# ------------------------------------------------
#  Figure C1: Actual and Estimated Slave Tax Rate
# ------------------------------------------------
library(ggpubr)
figC1 <- ggscatter(data, x = "actual_tax_over_value", y = "estimated_tax_over_value", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Actual Tax Rate", ylab = "Estimated Tax Rate")
figC1

# ------------------------------------
#  Figure D1: see file "Figure_D1.do"
# ------------------------------------

# -----------------------------------------------------------------------------------------------
#  Table G1: Cross-County Models: Enslaved Population and State Taxes per White Capita across Counties, 1860
# -----------------------------------------------------------------------------------------------
ctx <- read.csv("state_tax_counties.csv")

c0 <- lm(state_tax_pWHc ~ log(slaves_1860+1) + as.factor(statefips), data = ctx[which(ctx$malapportion == 1),])
se.c0 <- coeftest(c0, vcov=function(x) cluster.vcov(x, ctx[which(ctx$malapportion == 1),"statefips"]))[,2]

c1 <- lm(state_tax_pWHc ~ log(slaves_1860 +1) + ln_tot_pop + log(urban_pop_2500+1) + manu_gdp_pc_1860  + farm_value_pc_1860 + landgini_1860 + ln_area + as.factor(statefips), data = ctx[which(ctx$malapportion == 1),])
se.c1 <- coeftest(c1, vcov=function(x) cluster.vcov(x, ctx[which(ctx$malapportion == 1),"statefips"]))[,2]

c2 <- lm(state_tax_pWHc ~ log(slaves_1860 +1)*malapportion + as.factor(statefips), data = ctx)
se.c2 <- coeftest(c2, vcov=function(x) cluster.vcov(x, ctx$statefips))[,2]

c3 <- lm(state_tax_pWHc ~ log(slaves_1860 +1)*malapportion + ln_tot_pop + log(urban_pop_2500+1) + manu_gdp_pc_1860  + farm_value_pc_1860 + landgini_1860 + ln_area + as.factor(statefips), data = ctx)
se.c3 <- coeftest(c3, vcov=function(x) cluster.vcov(x, ctx$statefips))[,2]

stargazer(c0, c1, c2 , c3,
          se = list(se.c0, se.c1, se.c2 , se.c3), 
          keep= c("slaves_186", "malapportion"),
          covariate.labels = c("Enslaved population (log)", "Malapportionment", "Enslaved population (log) $times$ Malapportionment"),
          title = "Cross-County Models: Enslaved Population and State Taxes per White Capita across Counties, 1860",
          dep.var.labels = "State Taxes PWC",
          omit.stat = c("f", "ser", "adj.rsq"), omit = c("Constant"),
          add.lines = list(c("States", "7", "7", "14", "14"),
                           c("State FE", "Yes", "Yes", "Yes", "Yes"),
                           c("Additional covariates", "Yes", "Yes", "Yes", "Yes")) )



